1 //#**************************************************************
2 //# filename: ANCFBeamShear3D.cpp
3 //#
4 //# author: Karin & Astrid
5 //#
6 //# generated: Nov. 2010
7 //# description: 3D ANCF beam element with shear deformation
8 //# comments: formulation is based on the following paper:
9 //# "A 3D Shear Deformable Finite Element Based on the Absolute Nodal Coordinate Formulation"
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 "body3d.h"
25 #include "femathhelperfunctions.h"
26 #include "material.h"
27 #include "ANCFBeamShear3D.h"
28 #include "Node.h"
29 //#include "graphicsconstants.h"
30 #include "elementdataaccess.h"
31 //#include "solversettings_auto.h"
32
33
34
35
36 #pragma region ANCFBeamShear3DGeneric
37
CalculateElementLength() const38 double ANCFBeamShear3DGeneric::CalculateElementLength() const
39 {
40 //to obtain stress free reference configuration: same numerical integration as in axial deformation part of EvalF2()
41 int int_order_axial = (BeamFormulation == 4) ? 3 : order_axial; // must be set equally to integration order in axial deformation part of strain energy (EvalF2()) for each BeamFormulation, respectively!!!
42 Vector xT, wT;
43 GetIntegrationRule(xT, wT, int_order_axial);
44
45 double val = 0.;
46 for (int i1=1; i1<=xT.GetLen(); i1++) //i1-loop over all integration points
47 {
48 double x = xT(i1); //x in unit configuration [-1 , 1]
49 double dxidx = GetLx()*0.5; //determinant of the jacobian (from element transformation [-1 , 1] -> [-lx/2 , lx/2])
50 double xi = x*dxidx; //xi in scaled configuration [-lx/2 , lx/2]
51 Vector3D p0(xi, 0., 0.);
52
53 // integrate r,xi along xi in [-lx/2, +lx/2]
54 Vector sfx(NS(), 1); GetShapesx(sfx, p0);
55 Matrix Sfx(Dim(), NS()*Dim());
56 for (int i=1; i<=NS(); i++)
57 {
58 for (int j=1; j<=Dim(); j++)
59 {
60 Sfx(j, (i-1)*Dim() + j) = sfx(i);
61 }
62 }
63 Vector posx = Sfx*q0;
64 double norm_posx = posx.GetNorm();
65 val += norm_posx*wT(i1)*dxidx;
66 }
67
68 return val;
69 }
70
71
72 // default set function
SetANCFBeamShear3DGeneric(int materialnumi,const Vector3D & coli)73 void ANCFBeamShear3DGeneric::SetANCFBeamShear3DGeneric(int materialnumi, const Vector3D& coli)
74 {
75 //set reference configuration
76 int idx = 1;
77 for(int i=1; i <= nnodes; i++)
78 {
79 ANCFNodeS2S3_3D& n((ANCFNodeS2S3_3D&)GetNode(i));
80 q0(idx++) = n.Pos().X();
81 q0(idx++) = n.Pos().Y();
82 q0(idx++) = n.Pos().Z();
83 q0(idx++) = n.GetRefSlope2().X();
84 q0(idx++) = n.GetRefSlope2().Y();
85 q0(idx++) = n.GetRefSlope2().Z();
86 q0(idx++) = n.GetRefSlope3().X();
87 q0(idx++) = n.GetRefSlope3().Y();
88 q0(idx++) = n.GetRefSlope3().Z();
89 }
90
91 materialnum = materialnumi;
92
93 double lx = CalculateElementLength();
94 double ly, lz;
95 Beam3DProperties& bp((Beam3DProperties&)GetMaterial());
96 if (bp.IsMaterialOfBeamWithRectangularCrossSection())
97 {
98 ly = bp.GetBeamThicknessY();
99 lz = bp.GetBeamThicknessZ();
100 }
101 else
102 {
103 mbs->UO(UO_LVL_err) << "WARNING: Only rectangular cross section implemented for ANCFBeamShear3D at the moment!\n";
104 }
105 size.Set(lx,ly,lz);
106
107 xg.SetLen(SOS());
108 xg.SetAll(0.);
109
110 col = coli;
111 BuildDSMatrices();
112
113 penalty.SetLen(9);
114 penalty.SetAll(1.);
115
116 x_init.SetLen(2*SOS());
117 x_init.SetAll(0.);
118
119 SetElementName(GetElementSpec());
120 }
121
122 // deprecated set function!!!
SetANCFBeamShear3DGeneric(int materialnumi,const Vector3D & si,const Vector3D & coli)123 void ANCFBeamShear3DGeneric::SetANCFBeamShear3DGeneric(int materialnumi, const Vector3D& si, const Vector3D& coli)
124 {
125 size = si;
126
127 xg.SetLen(SOS());
128 xg.SetAll(0.);
129
130 //lx = size.X(); ly = size.Y(); lz = size.Z(); //lx is parameter to identify reference configuration
131 materialnum = materialnumi;
132 //mass = lx*ly*lz*GetMaterial().Density();
133 col = coli;
134 BuildDSMatrices();
135
136 penalty.SetLen(9);
137 penalty.SetAll(1.);
138
139 x_init.SetLen(2*SOS());
140 x_init.SetAll(0.);
141
142 SetElementName(GetElementSpec());
143 }
144
GetElementData(ElementDataContainer & edc)145 void ANCFBeamShear3DGeneric::GetElementData(ElementDataContainer& edc) //fill in all element data
146 {
147 Body3D::GetElementData(edc);
148 }
149
150 //user has changed data ==> write new data and parameters to element
SetElementData(ElementDataContainer & edc)151 int ANCFBeamShear3DGeneric::SetElementData(ElementDataContainer& edc) //fill in all element data
152 {
153 int rv = Body3D::SetElementData(edc);
154 return rv;
155 }
156
BuildDSMatrices()157 void ANCFBeamShear3DGeneric::BuildDSMatrices()
158 {
159 };
160
161
Gradu(const Vector3D & ploc,const Vector & u,Matrix3D & gradu) const162 void ANCFBeamShear3DGeneric::Gradu(const Vector3D& ploc, const Vector& u, Matrix3D& gradu) const
163 {
164 gradu.SetSize(3,3);
165 gradu.SetAll(0);
166
167 int dim = Dim();
168 int l;
169 for (int j = 1; j <= dim; j++)
170 {
171 for (int i = 1; i <= NS(); i++)
172 {
173 l = (i-1)*dim+j;
174 gradu(j,1) += GetSFx(i,ploc)*u(l);
175 gradu(j,2) += GetSFy(i,ploc)*u(l);
176 gradu(j,3) += GetSFz(i,ploc)*u(l);
177 }
178 }
179 }
180
181 //----------------
182 //Velocity
183 //----------------
184
185 //flagD = 0 f�r Berechnungen, flagD = 1 for Visualization
186
GetVel3D(const Vector3D & p_loc,int flagD) const187 Vector3D ANCFBeamShear3DGeneric::GetVel3D(const Vector3D& p_loc, int flagD) const
188 {
189 Vector3D p(0.,0.,0.);
190 for (int i = 1; i <= Dim(); i++) //i=1,2,3
191 {
192 for (int j = 1; j <= NS(); j++)
193 {
194 if(flagD==0)
195 p(i) += GetSF(j,p_loc)*XGP((j-1)*Dim()+i); //XGP=actual solution vector
196
197 else
198 p(i) += GetSF(j,p_loc)*XGPD((j-1)*Dim()+i);
199 }
200 }
201 return p;
202 };
203
204 //----------------
205 //Displacement
206 //----------------
207 //only for visualization
GetDisplacementD(const Vector3D & p_loc) const208 Vector3D ANCFBeamShear3DGeneric::GetDisplacementD(const Vector3D& p_loc) const
209 {
210 Vector3D p(0.,0.,0.);
211 for (int i = 1; i <= Dim(); i++)
212 {
213 for (int j = 1; j <= NS(); j++)
214 {
215 p(i) += GetSF(j,p_loc)*(XGD((j-1)*Dim()+i));
216 }
217 }
218 return p;
219 };
220
GetDisplacement(const Vector3D & p_loc) const221 Vector3D ANCFBeamShear3DGeneric::GetDisplacement(const Vector3D& p_loc) const
222 {
223 Vector3D p(0.,0.,0.);
224 for (int i = 1; i <= Dim(); i++)
225 {
226 for (int j = 1; j <= NS(); j++)
227 {
228 p(i) += GetSF(j,p_loc)*(XG((j-1)*Dim()+i));
229 }
230 }
231 return p;
232 };
233
234 //Vector3D ANCFBeamShear3DGeneric::GetDisplacement3D(const Vector3D& p_loc, const Vector& xg) const
235 //{
236 // Vector3D p(0.,0.,0.);
237 // for (int i = 1; i <= Dim(); i++)
238 // {
239 // for (int j = 1; j <= NS(); j++)
240 // {
241 // p(i) += GetSF(j,p_loc)*(xg((j-1)*Dim()+i));
242 // }
243 // }
244 // return p;
245 //};
246 //
247
GetIntDuDq(Matrix & dudq)248 void ANCFBeamShear3DGeneric::GetIntDuDq(Matrix& dudq)
249 {
250 //integrate du/dq over beam volume => (crossection area) * (integration over axis)
251 //u(q) = r(q) - r0 ==> du/dq = dr/dq
252
253 dudq.SetSize(NS()*Dim(), Dim());
254 dudq.SetAll(0);
255
256 ConstVector<4> w1, x1;//max. 4 integration points -> order 7
257 GetIntegrationRule(x1, w1, 3);
258 double fact = GetLx()*0.5*GetLy()*GetLz(); //crosssection area times determinant for line integral
259
260 for (int i1=1; i1<=x1.GetLen(); i1++)
261 {
262 double scaled_weight = w1(i1) * fact; //scale weight
263 double xi = 0.5*GetLx()*x1(i1);
264 Vector3D ploc(xi, 0., 0.);
265
266 for (int j=1; j<=NS(); j++)
267 {
268 double sf_value = GetSF(j, ploc); //value returned by i'th (positional) shape function at point xi
269 double block_j = Dim()*(j-1); //addresses (j-1)st block in dudq
270 for (int k=1; k<=Dim(); k++)
271 {
272 dudq(k+block_j, k) += scaled_weight*sf_value;
273 }
274 }
275 }
276 }
277
278 //----------------
279 //GetPos
280 //p_loc in [-Lx/2,Lx/2] x [-Ly/2,Ly/2] x [-Lz/2,Lz/2]
281 //----------------
282
GetPosD(const Vector3D & p_loc) const283 Vector3D ANCFBeamShear3DGeneric::GetPosD(const Vector3D& p_loc) const
284 {
285 return GetPos3D(p_loc, 1);
286 }
287
288 //Vector3D ANCFBeamShear3DGeneric::GetPosD(const Vector3D& p_loc) const
289 //// p_loc in [-Lx/2,Lx/2] x [-Ly/2,Ly/2] x [-Lz/2,Lz/2]
290 //{
291 // //mbs->UO() << "posy=" << GetPosy3D(p_loc, 1) << "\n";
292 // //mbs->UO() << "posz=" << GetPosz3D(p_loc, 1) << "\n";
293 // //mbs->UO() << "ploc=" << p_loc << "\n";
294 // Vector3D p0(p_loc.X(),0.,0.);
295 //
296 // return GetPos3D(p_loc, 1) + p_loc.Y()*GetPosy3D(p_loc, 1) + p_loc.Z()*GetPosz3D(p_loc, 1);
297 //}
298
299 //GetPos2D berechnet Vektor Positionsvektor r im deformierten Element
300 //GetPos2D computes r in deformed element
301 //xg = generalized coordinates (= our unknowns q)
GetPos3D(const Vector3D & p_loc,const Vector & xg) const302 Vector3D ANCFBeamShear3DGeneric::GetPos3D(const Vector3D& p_loc, const Vector& xg) const
303 {
304 Vector3D p(0.,0.,0.);
305 //TMStartTimer(24);
306 Vector sf(NS(), 1);
307 GetShapes(sf, p_loc);
308
309 for (int j = 1; j <= NS(); j++)
310 {
311 int offset = (j-1)*Dim();
312 for (int i = 1; i <= Dim(); i++)
313 {
314 p(i) += sf(j)*(xg(offset + i) + q0(offset + i)); //r,y=S,y*u+S,y*q_0
315 }
316 }
317 //TMStopTimer(24);
318 return p;
319 };
320
321 //GetPos2D von oben, aber falls im Aufruf als zweite Eingabe kein Vektor, sondern ein Integer steht,
322 //wird diese Fkt f�r Grafik aufgerufen; globale Variable XG wird f�r Rechnung verwendet anstatt Eingabe xg
323 //xg=XG wird einmal gesetzt, weshalb GetPos3D mit xg schneller ist, wohingegen GetPos3D mit XG jedes Mal aufs globale XG zugreift und das dauert lange
324
325 //same GetPos2D function as above, but with the second parameter flag D the function for the visualization is activated
GetPos3D(const Vector3D & p_loc,int flagD) const326 Vector3D ANCFBeamShear3DGeneric::GetPos3D(const Vector3D& p_loc, int flagD) const
327 {
328 Vector3D p(0.,0.,0.);
329 Vector sf(NS(), 1);
330 GetShapes(sf, p_loc);
331
332 if (!flagD)
333 {
334 for (int j = 1; j <= NS(); j++)
335 {
336 int offset = (j-1)*Dim();
337 for (int i = 1; i <= Dim(); i++)
338 {
339 p(i) += sf(j)*(XG(offset + i) + q0(offset + i));
340 }
341 }
342 }
343 else
344 {
345 for (int j = 1; j <= NS(); j++)
346 {
347 int offset = (j-1)*Dim();
348 for (int i = 1; i <= Dim(); i++)
349 {
350 p(i) += sf(j)*(XGD(offset + i) + q0(offset + i));
351 }
352 }
353 }
354 return p;
355 };
356
357 //Ableitungen von GetPos entsprechen unseren Richtungsvektoren
358 //GetPosx2D (Ableitung in Richtung der Balkenachse: dr/dx=r,x)
359 //GetPosx2D is the derivative with respect to the direction of the beam centerline
GetPosx3D(const Vector3D & p_loc,const Vector & xg) const360 Vector3D ANCFBeamShear3DGeneric::GetPosx3D(const Vector3D& p_loc, const Vector& xg) const
361 {
362 Vector3D p(0.,0.,0.);
363 Vector sfx(NS(), 1);
364 GetShapesx(sfx, p_loc);
365
366 for (int j = 1; j <= NS(); j++)
367 {
368 int offset = (j-1)*Dim();
369 for (int i = 1; i <= Dim(); i++)
370 {
371 p(i) += sfx(j)*(xg(offset + i) + q0(offset + i)); //r,y=S,y*u+S,y*q_0
372 }
373 }
374 return p;
375 }
376
GetPosx3D(const Vector3D & p_loc,int flagD) const377 Vector3D ANCFBeamShear3DGeneric::GetPosx3D(const Vector3D& p_loc, int flagD) const
378 {
379 Vector3D p(0.,0.,0.);
380 Vector sfx(NS(), 1);
381 GetShapesx(sfx, p_loc);
382
383 if (!flagD)
384 {
385 for (int j = 1; j <= NS(); j++)
386 {
387 int offset = (j-1)*Dim();
388 for (int i = 1; i <= Dim(); i++)
389 {
390 p(i) += sfx(j)*(XG(offset + i) + q0(offset + i));
391 }
392 }
393 }
394 else
395 {
396 for (int j = 1; j <= NS(); j++)
397 {
398 int offset = (j-1)*Dim();
399 for (int i = 1; i <= Dim(); i++)
400 {
401 p(i) += sfx(j)*(XGD(offset + i) + q0(offset + i));
402 }
403 }
404 }
405 return p;
406 }
407
408
GetdPosdx(const Vector3D & ploc,Vector3D & dpdx)409 void ANCFBeamShear3DGeneric::GetdPosdx(const Vector3D& ploc, Vector3D& dpdx)
410 {
411 dpdx = GetPosx3D(ploc);
412 };
413
414
415 //GetPosy2D (Ableitung in Richtung des Querschnitts: dr/dy=r,y)
416 //GetPosy2D is the derivative with respect to the direction of the cross section
GetPosy3D(const Vector3D & p_loc,const Vector & xg) const417 Vector3D ANCFBeamShear3DGeneric::GetPosy3D(const Vector3D& p_loc, const Vector& xg) const
418 {
419 Vector3D p(0.,0.,0.);
420 Vector sfy(NS(), 1);
421 GetShapesy(sfy, p_loc);
422
423 for (int j = 1; j <= NS(); j++)
424 {
425 int offset = (j-1)*Dim();
426 for (int i = 1; i <= Dim(); i++)
427 {
428 p(i) += sfy(j)*(xg(offset + i) + q0(offset + i)); //r,y=S,y*u+S,y*q_0
429 }
430 }
431 return p;
432 }
433
GetPosy3D(const Vector3D & p_loc,int flagD) const434 Vector3D ANCFBeamShear3DGeneric::GetPosy3D(const Vector3D& p_loc, int flagD) const
435 {
436 Vector3D p(0.,0.,0.);
437 Vector sfy(NS(), 1);
438 GetShapesy(sfy, p_loc);
439 if (!flagD)
440 {
441 for (int j = 1; j <= NS(); j++)
442 {
443 int offset = (j-1)*Dim();
444 for (int i = 1; i <= Dim(); i++)
445 {
446 p(i) += sfy(j)*(XG(offset + i) + q0(offset + i));
447 }
448 }
449 }
450 else
451 {
452 for (int j = 1; j <= NS(); j++)
453 {
454 int offset = (j-1)*Dim();
455 for (int i = 1; i <= Dim(); i++)
456 {
457 p(i) += sfy(j)*(XGD(offset + i) + q0(offset + i));
458 }
459 }
460 }
461 return p;
462 }
463
GetPosyP3D(const Vector3D & p_loc,int flagD) const464 Vector3D ANCFBeamShear3DGeneric::GetPosyP3D(const Vector3D& p_loc, int flagD) const
465 {
466 Vector3D p(0.,0.,0.);
467 Vector sfy(NS(), 1);
468 GetShapesy(sfy, p_loc);
469 if (!flagD)
470 {
471 for (int j = 1; j <= NS(); j++)
472 {
473 int offset = (j-1)*Dim();
474 for (int i = 1; i <= Dim(); i++)
475 {
476 p(i) += sfy(j)*XGP(offset + i);
477 }
478 }
479 }
480 else
481 {
482 for (int j = 1; j <= NS(); j++)
483 {
484 int offset = (j-1)*Dim();
485 for (int i = 1; i <= Dim(); i++)
486 {
487 p(i) += sfy(j)*XGPD(offset + i);
488 }
489 }
490 }
491 return p;
492 }
493
GetPosz3D(const Vector3D & p_loc,const Vector & xg) const494 Vector3D ANCFBeamShear3DGeneric::GetPosz3D(const Vector3D& p_loc, const Vector& xg) const
495 {
496 Vector3D p(0.,0.,0);
497 Vector sfz(NS(), 1);
498 GetShapesz(sfz, p_loc);
499
500 for (int j = 1; j <= NS(); j++)
501 {
502 int offset = (j-1)*Dim();
503 for (int i = 1; i <= Dim(); i++)
504 {
505 p(i) += sfz(j)*(xg(offset + i) + q0(offset + i)); //r,y=S,y*u+S,y*q_0
506 }
507 }
508 return p;
509 }
510
GetPosz3D(const Vector3D & p_loc,int flagD) const511 Vector3D ANCFBeamShear3DGeneric::GetPosz3D(const Vector3D& p_loc, int flagD) const
512 {
513 Vector3D p(0.,0.,0.);
514 Vector sfz(NS(), 1);
515 GetShapesz(sfz, p_loc);
516
517 if (!flagD)
518 {
519 for (int j = 1; j <= NS(); j++)
520 {
521 int offset = (j-1)*Dim();
522 for (int i = 1; i <= Dim(); i++)
523 {
524 p(i) += sfz(j)*(XG(offset + i) + q0(offset + i));
525 }
526 }
527 }
528 else
529 {
530 for (int j = 1; j <= NS(); j++)
531 {
532 int offset = (j-1)*Dim();
533 for (int i = 1; i <= Dim(); i++)
534 {
535 p(i) += sfz(j)*(XGD(offset + i) + q0(offset + i));
536 }
537 }
538 }
539 return p;
540 }
541
GetPoszP3D(const Vector3D & p_loc,int flagD) const542 Vector3D ANCFBeamShear3DGeneric::GetPoszP3D(const Vector3D& p_loc, int flagD) const
543 {
544 Vector3D p(0.,0.,0.);
545 Vector sfz(NS(), 1);
546 GetShapesz(sfz, p_loc);
547
548 if (!flagD)
549 {
550 for (int j = 1; j <= NS(); j++)
551 {
552 int offset = (j-1)*Dim();
553 for (int i = 1; i <= Dim(); i++)
554 {
555 p(i) += sfz(j)*XGP(offset + i);
556 }
557 }
558 }
559 else
560 {
561 for (int j = 1; j <= NS(); j++)
562 {
563 int offset = (j-1)*Dim();
564 for (int i = 1; i <= Dim(); i++)
565 {
566 p(i) += sfz(j)*XGPD(offset + i);
567 }
568 }
569 }
570 return p;
571 }
572
573
574 //GetPosxy2D: d2r/dydx=r,xy
575 //is used in DeltaThetax
GetPosxy3D(const Vector3D & p_loc,const Vector & xg) const576 Vector3D ANCFBeamShear3DGeneric::GetPosxy3D(const Vector3D& p_loc, const Vector& xg) const
577 {
578 Vector3D p(0.,0.,0.);
579 Vector sfxy(NS(), 1);
580 GetShapesxy(sfxy, p_loc);
581
582 for (int j = 1; j <= NS(); j++)
583 {
584 int offset = (j-1)*Dim();
585 for (int i = 1; i <= Dim(); i++)
586 {
587 p(i) += sfxy(j)*(xg(offset + i) + q0(offset + i)); //r,y=S,y*u+S,y*q_0
588 }
589 }
590 return p;
591 }
592
GetPosxy3D(const Vector3D & p_loc,int flagD) const593 Vector3D ANCFBeamShear3DGeneric::GetPosxy3D(const Vector3D& p_loc, int flagD) const
594 {
595 Vector3D p(0.,0.,0.);
596 Vector sfxy(NS(), 1);
597 GetShapesxy(sfxy, p_loc);
598
599 if (!flagD)
600 {
601 for (int j = 1; j <= NS(); j++)
602 {
603 int offset = (j-1)*Dim();
604 for (int i = 1; i <= Dim(); i++)
605 {
606 p(i) += sfxy(j)*(XG(offset + i) + q0(offset + i));
607 }
608 }
609 }
610 else
611 {
612 for (int j = 1; j <= NS(); j++)
613 {
614 int offset = (j-1)*Dim();
615 for (int i = 1; i <= Dim(); i++)
616 {
617 p(i) += sfxy(j)*(XGD(offset + i) + q0(offset + i));
618 }
619 }
620 }
621 return p;
622 }
623
GetPosxz3D(const Vector3D & p_loc,const Vector & xg) const624 Vector3D ANCFBeamShear3DGeneric::GetPosxz3D(const Vector3D& p_loc, const Vector& xg) const
625 {
626 Vector3D p(0.,0.,0.);
627 Vector sfxz(NS(), 1);
628 GetShapesxz(sfxz, p_loc);
629
630 for (int j = 1; j <= NS(); j++)
631 {
632 int offset = (j-1)*Dim();
633 for (int i = 1; i <= Dim(); i++)
634 {
635 p(i) += sfxz(j)*(xg(offset + i) + q0(offset + i)); //r,y=S,y*u+S,y*q_0
636 }
637 }
638 return p;
639 }
640
GetPosxz3D(const Vector3D & p_loc,int flagD) const641 Vector3D ANCFBeamShear3DGeneric::GetPosxz3D(const Vector3D& p_loc, int flagD) const
642 {
643 Vector3D p(0.,0.,0.);
644 Vector sfxz(NS(), 1);
645 GetShapesxz(sfxz, p_loc);
646
647 if (!flagD)
648 {
649 for (int j = 1; j <= NS(); j++)
650 {
651 int offset = (j-1)*Dim();
652 for (int i = 1; i <= Dim(); i++)
653 {
654 p(i) += sfxz(j)*(XG(offset + i) + q0(offset + i));
655 }
656 }
657 }
658 else
659 {
660 for (int j = 1; j <= NS(); j++)
661 {
662 int offset = (j-1)*Dim();
663 for (int i = 1; i <= Dim(); i++)
664 {
665 p(i) += sfxz(j)*(XGD(offset + i) + q0(offset + i));
666 }
667 }
668 }
669 return p;
670 }
671
672 //computes position vector in the reference element: r_0
673 //p_loc in [-Lx/2,Lx/2] x [-Ly/2,Ly/2] x [-Lz/2,Lz/2]
GetInitPos3D(const Vector3D & p_loc) const674 Vector3D ANCFBeamShear3DGeneric::GetInitPos3D(const Vector3D& p_loc) const
675 {
676 Vector3D p(0.,0.,0.);
677 Vector sf(NS(), 1);
678 GetShapes(sf, p_loc);
679 for (int j = 1; j <= NS(); j++)
680 {
681 int offset = (j-1)*Dim();
682 for (int i = 1; i <= Dim(); i++)
683 {
684 p(i) += sf(j)*q0(offset + i);
685 }
686 }
687 return p;
688 };
689
GetInitPosx3D(const Vector3D & p_loc) const690 Vector3D ANCFBeamShear3DGeneric::GetInitPosx3D(const Vector3D& p_loc) const
691 {
692 Vector3D p(0.,0.,0.);
693 Vector sf(NS(), 1);
694 GetShapesx(sf, p_loc);
695 for (int j = 1; j <= NS(); j++)
696 {
697 int offset = (j-1)*Dim();
698 for (int i = 1; i <= Dim(); i++)
699 {
700 p(i) += sf(j)*q0(offset + i);
701 }
702 }
703 return p;
704 };
705
GetInitPosy3D(const Vector3D & p_loc) const706 Vector3D ANCFBeamShear3DGeneric::GetInitPosy3D(const Vector3D& p_loc) const
707 {
708 Vector3D p(0.,0.,0.);
709 Vector sf(NS(), 1);
710 GetShapesy(sf, p_loc);
711 for (int j = 1; j <= NS(); j++)
712 {
713 int offset = (j-1)*Dim();
714 for (int i = 1; i <= Dim(); i++)
715 {
716 p(i) += sf(j)*q0(offset + i);
717 }
718 }
719 return p;
720 }
721
GetInitPosz3D(const Vector3D & p_loc) const722 Vector3D ANCFBeamShear3DGeneric::GetInitPosz3D(const Vector3D& p_loc) const
723 {
724 Vector3D p(0.,0.,0.);
725 Vector sf(NS(), 1);
726 GetShapesz(sf, p_loc);
727 for (int j = 1; j <= NS(); j++)
728 {
729 int offset = (j-1)*Dim();
730 for (int i = 1; i <= Dim(); i++)
731 {
732 p(i) += sf(j)*q0(offset + i);
733 }
734 }
735 return p;
736 }
737
GetInitPosxy3D(const Vector3D & p_loc) const738 Vector3D ANCFBeamShear3DGeneric::GetInitPosxy3D(const Vector3D& p_loc) const
739 {
740 Vector3D p(0.,0.,0.);
741 Vector sf(NS(), 1);
742 GetShapesxy(sf, p_loc);
743 for (int j = 1; j <= NS(); j++)
744 {
745 int offset = (j-1)*Dim();
746 for (int i = 1; i <= Dim(); i++)
747 {
748 p(i) += sf(j)*q0(offset + i);
749 }
750 }
751 return p;
752 }
753
GetInitPosxz3D(const Vector3D & p_loc) const754 Vector3D ANCFBeamShear3DGeneric::GetInitPosxz3D(const Vector3D& p_loc) const
755 {
756 Vector3D p(0.,0.,0.);
757 Vector sf(NS(), 1);
758 GetShapesxz(sf, p_loc);
759 for (int j = 1; j <= NS(); j++)
760 {
761 int offset = (j-1)*Dim();
762 for (int i = 1; i <= Dim(); i++)
763 {
764 p(i) += sf(j)*q0(offset + i);
765 }
766 }
767 return p;
768 }
769
770 //for visualization
GetPos3D0D(const Vector3D & p_loc) const771 Vector3D ANCFBeamShear3DGeneric::GetPos3D0D(const Vector3D& p_loc) const
772 {
773 Vector3D plocscaled;
774 plocscaled(1)=p_loc(1)*GetLx()/2.;
775 plocscaled(2)=p_loc(2)*GetLy()/2.;
776 plocscaled(3)=p_loc(3)*GetLz()/2.;
777 return GetPos3D(plocscaled,1);
778 };
779
GetPos3D0D(const Vector3D & p_loc,double defscale) const780 Vector3D ANCFBeamShear3DGeneric::GetPos3D0D(const Vector3D& p_loc, double defscale) const
781 {
782 Vector3D plocscaled;
783 plocscaled(1)=p_loc(1)*GetLx()/2.;
784 plocscaled(2)=p_loc(2)*GetLy()/2.;
785 plocscaled(3)=p_loc(3)*GetLz()/2.;
786 return GetInitPos3D(plocscaled) + defscale*GetDisplacementD(plocscaled);
787
788 GetNodePos(1);
789 };
790
791 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
GetNodeLocPos(int i) const792 Vector3D ANCFBeamShear3DGeneric::GetNodeLocPos(int i) const //returns position of i-th node
793 {
794 switch(i)
795 {
796 case 1: return Vector3D(-GetLx()*.5, 0, 0);
797 case 2: return Vector3D(+GetLx()*.5, 0, 0);
798 case 3: assert(nnodes >= 3); return Vector3D(0.);
799 }
800 assert(0);
801 return Vector3D(0.);
802 }
803
804
GetNodePos(int i) const805 Vector3D ANCFBeamShear3DGeneric::GetNodePos(int i) const //returns position of i-th node
806 {
807 switch(i)
808 {
809 case 1: return GetPos(Vector3D(-GetLx()*.5, 0, 0));
810 case 2: return GetPos(Vector3D(+GetLx()*.5, 0, 0));
811 case 3: assert(nnodes >= 3); return GetPos(Vector3D(0.));
812 }
813 assert(0);
814 return Vector3D(0.);
815 }
816
GetNodePosD(int i) const817 Vector3D ANCFBeamShear3DGeneric::GetNodePosD(int i) const //returns position of i-th node (draw mode)
818 {
819 switch(i)
820 {
821 case 1: return GetPosD(Vector3D(-GetLx()*.5, 0, 0));
822 case 2: return GetPosD(Vector3D(+GetLx()*.5, 0, 0));
823 case 3: assert(nnodes >= 3); return GetPosD(Vector3D(0.));
824 }
825 assert(0);
826 return Vector3D(0.);
827 }
828 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
829
GetDOFDirD(int idof) const830 Vector3D ANCFBeamShear3DGeneric::GetDOFDirD(int idof) const //returns direction of action of i-th DOF
831 {
832 if ((idof >= 1 && idof <= 3) ||
833 (idof >= 10 && idof <= 12) ||
834 (idof >= 19 && idof <= 21))
835 {
836 Vector3D v(0.);
837 int component = (idof-1)%3+1;
838 v(component) = 1;
839 return v;
840 }
841 //for other components, the direction can not be easily identified and is therefore set to zero (default)
842
843 return Vector3D(0.,0.,0.);
844 }
845
GetDOFPosD(int idof) const846 Vector3D ANCFBeamShear3DGeneric::GetDOFPosD(int idof) const //returns position of i-th DOF
847 {
848 if (idof <= 9) //left node: 1,...,9
849 {
850 return GetPos3D(Vector3D(-0.5*GetLx(),0.,0.),1);//1 f�r flagD
851 }
852 else if (idof <= 18) //right node: 10,...,18
853 {
854 return GetPos3D(Vector3D(0.5*GetLx(),0.,0.),1);
855 }
856 else if (idof <= 27) //middle node: 19,...,27
857 {
858 return GetPos3D(Vector3D(0.,0.,0.),1);//middle node: (0,0)
859 }
860 else
861 {
862 return Vector3D(0, 0, 0);
863 }
864 }
865
866
GetRotMatrix(const Vector3D & p0,int on_reference_element,int flagD) const867 Matrix3D ANCFBeamShear3DGeneric::GetRotMatrix(const Vector3D& p0, int on_reference_element, int flagD) const //p0(1) in [-lx/2,lx/2]
868 {
869 //ROTATION TENSOR:
870 //A = [ e1 | e2 | e3 ]
871
872 //e1 = e1^{\bar}/|e1^{\bar}| with e1^{\bar} = drdy x drdz
873 //e3 = e3^{\bar}/|e3^{\bar}| with e3^{\bar} = drdz
874 //e2 = e2^{\bar}/|e2^{\bar}| with e2^{\bar} = drdz x (drdy x drdz) = e3^{\bar} x e1^{\bar}
875
876 //Attention: Getei(i) returns ei^{\bar}
877
878 Vector3D drdy, drdz;
879
880 if (on_reference_element)
881 {
882 drdy = GetInitPosy3D(p0);
883 drdz = GetInitPosz3D(p0);
884 }
885 else
886 {
887 drdy = GetPosy3D(p0, flagD);
888 drdz = GetPosz3D(p0, flagD);
889 }
890
891 TArray<Vector3D> e(Dim());
892 for (int i=1; i<=Dim(); i++)
893 {
894 e(i) = Getei(i, drdy.X(), drdy.Y(), drdy.Z(), drdz.X(), drdz.Y(), drdz.Z());
895 e(i).Normalize();
896 }
897
898 Matrix3D A(e(1).X(),e(2).X(),e(3).X(),
899 e(1).Y(),e(2).Y(),e(3).Y(),
900 e(1).Z(),e(2).Z(),e(3).Z());
901 return A;
902 }
903
GetRotMatrixP(const Vector3D & p0,int flagD) const904 Matrix3D ANCFBeamShear3DGeneric::GetRotMatrixP(const Vector3D& p0, int flagD) const
905 {
906 //Time derivative of rotation tensor:
907 //AP = [ e1P | e2P | e3P ]
908
909 //e1 = e1^{\bar}/|e1^{\bar}| with eP1^{\bar} = drdyP x drdz + drdy x drdzP
910 //e3 = e3^{\bar}/|e3^{\bar}| with eP3^{\bar} = drdzP
911 //e2 = e2^{\bar}/|e2^{\bar}| with eP2^{\bar} = drdzP x (drdy x drdz) + drdz x (drdyP x drdz + drdy x drdzP) = eP3^{\bar} x e1^{\bar} + e3^{\bar} x eP1^{\bar}
912
913 //Attention: Getei(i) returns ei^{\bar}, and GetePi(i) returns ePi^{\bar}
914
915 Vector3D drdy, drdz, drPdy, drPdz;
916
917 drdy = GetPosy3D(p0, flagD);
918 drdz = GetPosz3D(p0, flagD);
919 drPdy = GetPosyP3D(p0, flagD);
920 drPdz = GetPoszP3D(p0, flagD);
921
922 Matrix3D AP;
923 for (int i=1; i<=3; i++)
924 {
925 Vector3D ei = Getei(i, drdy.X(), drdy.Y(), drdy.Z(), drdz.X(), drdz.Y(), drdz.Z());
926 double ei_norm = ei.Norm();
927 Vector3D eiP = GeteiP(i, drdy.X(), drdy.Y(), drdy.Z(), drdz.X(), drdz.Y(), drdz.Z(), drPdy.X(), drPdy.Y(), drPdy.Z(), drPdz.X(), drPdz.Y(), drPdz.Z());
928 double ei_eiP = ei*eiP;
929 eiP *= (1./ei_norm);
930 eiP += (-ei_eiP/pow(ei_norm,3))*ei;
931 for (int j=1; j<=3; j++)
932 {
933 AP(j,i) = eiP(j);
934 }
935 }
936
937 return AP;
938 }
939
GetTwistAndCurvature(const Vector3D & ploc,int on_reference_element,int flagD) const940 Vector3D ANCFBeamShear3DGeneric::GetTwistAndCurvature(const Vector3D& ploc, int on_reference_element, int flagD) const
941 // ploc.X() in [ -Lx/2 , Lx/2 ]
942 // calculate vector of twist and curvature k = 1/2 sum_i e_i x e_i'
943 // or, for the reference element (see flag) k_0 = 1/2 sum_i e0_i x e0_i'
944 // NOTE: both vectors are returned with their components w.r.t. the local frame basis (e_i)!
945 {
946 if(IsStraightBeamInReferenceConfiguration() && on_reference_element) // for speedup
947 {
948 return Vector3D(); //std constructor initializes with zero
949 }
950
951 Vector3D drdy, drdz, ddrdxdy, ddrdxdz;
952 if (on_reference_element)
953 {
954 drdy = GetInitPosy3D(ploc);
955 drdz = GetInitPosz3D(ploc);
956 ddrdxdy = GetInitPosxy3D(ploc);
957 ddrdxdz = GetInitPosxz3D(ploc);
958 }
959 else
960 {
961 drdy = GetPosy3D(ploc, flagD);
962 drdz = GetPosz3D(ploc, flagD);
963 ddrdxdy = GetPosxy3D(ploc, flagD);
964 ddrdxdz = GetPosxz3D(ploc, flagD);
965 }
966
967 Vector3D ei; // returns ei^{\bar}, ei = ei^{\bar} / |ei^{\bar}|
968 //e1 = e1^{\bar}/|e1^{\bar}| with e1^{\bar} = drdy x drdz
969 //e3 = e3^{\bar}/|e3^{\bar}| with e3^{\bar} = drdz
970 //e2 = e2^{\bar}/|e2^{\bar}| with e2^{\bar} = drdz x (drdy x drdz) = e3^{\bar} x e1^{\bar}
971 Vector3D deidxi; // returns d ei^{\bar} / d xi
972
973 Vector3D k;
974 for (int i=1; i<=Dim(); i++)
975 {
976 ei = Getei(i, drdy.X(), drdy.Y(), drdy.Z(), drdz.X(), drdz.Y(), drdz.Z());
977 deidxi = Getdeidxi(i, drdy.X(), drdy.Y(), drdy.Z(), drdz.X(), drdz.Y(), drdz.Z(), ddrdxdy.X(), ddrdxdy.Y(), ddrdxdy.Z(), ddrdxdz.X(), ddrdxdz.Y(), ddrdxdz.Z());
978
979 k += ( 0.5 / ei.Norm2() ) * ei.Cross(deidxi); //k = 1/2 sum_i e_i x e_i' -> simplified written using e_i^{\bar}-Vectors (see the paper)
980 }
981
982 Matrix3D AT = GetRotMatrix(ploc, on_reference_element, flagD).GetTp(); //in order to represent the vector k in the local frame basis
983
984 return AT * k;
985 }
986
GetAngularVel(const Vector3D & p_loc) const987 Vector3D ANCFBeamShear3DGeneric::GetAngularVel(const Vector3D& p_loc) const
988 {
989 Matrix retmat;
990 Vector xp(SOS());
991 for (int j=1; j<=SOS(); j++)
992 {
993 xp(j)=XGP(j);
994 }
995 Matrix3D A = GetRotMatrix(p_loc);
996 Matrix3D omega_skew;
997 for (int i=1; i<=Dim(); i++)
998 {
999 Vector3D ei(A(1,i),A(2,i),A(3,i));
1000 GetdRotvdqT(ei, p_loc, retmat, 1);
1001 Vector wi(3);
1002 Mult(xp, retmat, wi);
1003 for (int j=1; j<=Dim(); j++)
1004 {
1005 omega_skew(i,j) = wi(j);
1006 }
1007 }
1008 //Matrix3D omega_skew = GetRotMatrixP(p_loc)*GetRotMatrix(p_loc).GetTp();
1009 //omega_skew * A = RotMatrixP
1010
1011 return Vector3D(omega_skew(2,3),omega_skew(1,3),-omega_skew(1,2));
1012 }
1013
GetdRotvdqT(const Vector3D & vloc,const Vector3D & ploc,Matrix & drotvdqT,int use_transposed_rotmatrix) const1014 void ANCFBeamShear3DGeneric::GetdRotvdqT(const Vector3D& vloc, const Vector3D& ploc, Matrix& drotvdqT, int use_transposed_rotmatrix) const
1015 {
1016 // rot = (rot11 rot12 rot13
1017 // rot21 rot22 rot23
1018 // rot31 rot32 rot33)
1019 // Basis vectors are columns of rotation tensor
1020 // rotT.v = (rot11*v1 + rot21*v2 + rot31*v3, rot12*v1 + rot22*v2 + rot32*v3, rot13*v1 + rot23*v2 + rot33*v3)
1021 // p = (r,y , r,z) = ( r1dy , r2dy , r3dy , r1dz , r2dz , r3dz )
1022 // drotTv/dq = drotTv/dp . dp/dq
1023
1024 double xi = ploc.X(); // ploc.X() in [ -Lx/2 , Lx/2 ]
1025
1026 //rotation tensor only depends on r,y and r,z
1027 Vector3D drdy = GetPosy3D(xi);
1028 Vector3D drdz = GetPosz3D(xi);
1029
1030 Vector3D e1(Getei(1, drdy.X(), drdy.Y(), drdy.Z(), drdz.X(), drdz.Y(), drdz.Z()));
1031 Vector3D e2(Getei(2, drdy.X(), drdy.Y(), drdy.Z(), drdz.X(), drdz.Y(), drdz.Z()));
1032 Vector3D e3(Getei(3, drdy.X(), drdy.Y(), drdy.Z(), drdz.X(), drdz.Y(), drdz.Z()));
1033 TArray<Vector3D> de1dp(GetdeidDiffVar(1, drdy.X(), drdy.Y(), drdy.Z(), drdz.X(), drdz.Y(), drdz.Z()));
1034 TArray<Vector3D> de2dp(GetdeidDiffVar(2, drdy.X(), drdy.Y(), drdy.Z(), drdz.X(), drdz.Y(), drdz.Z()));
1035 TArray<Vector3D> de3dp(GetdeidDiffVar(3, drdy.X(), drdy.Y(), drdy.Z(), drdz.X(), drdz.Y(), drdz.Z()));
1036
1037 // calculate derivatives of normalized ei
1038 double e1norm = e1.Norm();
1039 double e2norm = e2.Norm();
1040 double e3norm = e3.Norm();
1041 for (int j=1; j<=de1dp.Length(); j++)
1042 {
1043 de1dp(j) = de1dp(j)*(1./e1norm) - e1*((e1*de1dp(j))/Cub(e1norm));
1044 de2dp(j) = de2dp(j)*(1./e2norm) - e2*((e2*de2dp(j))/Cub(e2norm));
1045 de3dp(j) = de3dp(j)*(1./e3norm) - e3*((e3*de3dp(j))/Cub(e3norm));
1046 }
1047
1048 Matrix drotvdp(Dim(), 6);//3x6
1049 if (use_transposed_rotmatrix) //use transposed A: A^T*v
1050 {
1051 for (int h=1; h<=6; h++)//h...DiffVar
1052 {
1053 drotvdp(1,h) = de1dp(h)*vloc;
1054 drotvdp(2,h) = de2dp(h)*vloc;
1055 drotvdp(3,h) = de3dp(h)*vloc;
1056 }
1057 }
1058 else //dont use transposed A: A*v
1059 {
1060 for (int h=1; h<=6; h++) //h...DiffVar
1061 {
1062 drotvdp(1,h) = de1dp(h).X()*vloc.X() + de2dp(h).X()*vloc.Y() + de3dp(h).X()*vloc.Z();
1063 drotvdp(2,h) = de1dp(h).Y()*vloc.X() + de2dp(h).Y()*vloc.Y() + de3dp(h).Y()*vloc.Z();
1064 drotvdp(3,h) = de1dp(h).Z()*vloc.X() + de2dp(h).Z()*vloc.Y() + de3dp(h).Z()*vloc.Z();
1065 }
1066 }
1067 //h = 1,2,3 -> SFy
1068 //h = 4,5,6 -> SFz
1069 drotvdqT.SetSize(SOS(),Dim());
1070 drotvdqT.SetAll(0.);
1071 Vector sfy(NS(), 1); GetShapesy(sfy, ploc);
1072 Vector sfz(NS(), 1); GetShapesz(sfz, ploc);
1073
1074 for (int k=1; k<=Dim(); k++)
1075 {
1076 for (int i=1; i<=3; i++)
1077 for (int j=1; j<=NS(); j++)
1078 drotvdqT((j-1)*Dim()+i, k) = drotvdp(k,i)*sfy(j) + drotvdp(k,Dim()+i)*sfz(j);
1079 }
1080 }
1081
GetdRotdqT(const Vector3D & ploc,Matrix & d)1082 void ANCFBeamShear3DGeneric::GetdRotdqT(const Vector3D& ploc, Matrix& d)
1083 {
1084 //ploc in [-lx/2,lx/2] x [-ly/2,ly/2] x [-lz/2,lz/2]
1085 ////////////////////////////////////////////////////////////
1086 //KN&PG's approach \cite{hodges06}
1087 //let e_i denote the i-th base vector of the actual system (i-th column of the rotation matrix)
1088 //\delta \theta = \sum_{i=1}^3 e_i (\delta e_j . e_k) where j = i%3+1, k = j%3+1
1089 //(\grad_q \theta^T)_{l,m} = (e_i)_m (e_k)_n (partial (e_j)_n/ \partial q_l\)
1090 //note: for this method the call of routine GetdRotvdqT has to sum up (e_j)_n v_n instead of (e_j)_n v_j, i.e. use_transposed_rotmatrix = 1
1091 ////////////////////////////////////////////////////////////
1092
1093 d.SetSize(SOS(),Dim());
1094 d.FillWithZeros();
1095
1096 double xi = ploc.X(); // ploc.X() in [ -Lx/2 , Lx/2 ]
1097
1098 //rotation tensor only depends on r,y and r,z
1099 Vector3D drdy = GetPosy3D(xi);
1100 Vector3D drdz = GetPosz3D(xi);
1101
1102 //-----------------------------------------
1103 // new version
1104 //-----------------------------------------
1105 Vector3D ei;
1106 Vector3D ej;
1107 Vector3D ek;
1108 TArray<Vector3D> dejdp(6);
1109 TArray<Vector3D> dekdp(6);
1110 TArray<Vector3D> deidp(6);
1111 ConstVector<27> helpy;
1112 helpy.SetAll(0.);
1113
1114
1115 bool define_theta_by_a_transposed = false;
1116 bool derivatives_wrt_rotated_vars = false;
1117
1118 for (int i=1; i<=Dim(); i++) //i=1,2,3
1119 {
1120 int j = i%3 + 1; //i=1->j=2->k=3, i=2->j=3->k=1, i=3->j=1->k=2
1121 int k = j%3 + 1;
1122
1123 // e_i = ebar_i / |ebar_i|
1124 ei = Getei(i, drdy.X(), drdy.Y(), drdy.Z(), drdz.X(), drdz.Y(), drdz.Z());
1125 ej = Getei(j, drdy.X(), drdy.Y(), drdy.Z(), drdz.X(), drdz.Y(), drdz.Z());
1126 ek = Getei(k, drdy.X(), drdy.Y(), drdy.Z(), drdz.X(), drdz.Y(), drdz.Z());
1127
1128 if (define_theta_by_a_transposed)
1129 {
1130 double memo;
1131 memo = ei(j);
1132 ei(j) = ej(i);
1133 ej(i) = memo;
1134 memo = ej(k);
1135 ej(k) = ek(j);
1136 ek(j) = memo;
1137 memo = ek(i);
1138 ek(i) = ei(k);
1139 ei(k) = memo;
1140 }
1141
1142 ei.Normalize();
1143 ek.Normalize();
1144 double ej_norm = ej.Norm();
1145 double ej_norm3 = Cub(ej_norm);
1146 ej.Normalize();
1147
1148 // GetdeidDiffVar(j,...) returns debar_j_dp
1149 deidp = GetdeidDiffVar(i, drdy.X(), drdy.Y(), drdy.Z(), drdz.X(), drdz.Y(), drdz.Z());
1150 dejdp = GetdeidDiffVar(j, drdy.X(), drdy.Y(), drdy.Z(), drdz.X(), drdz.Y(), drdz.Z());
1151 dekdp = GetdeidDiffVar(k, drdy.X(), drdy.Y(), drdy.Z(), drdz.X(), drdz.Y(), drdz.Z());
1152
1153 // calculate all n=2*3=6 partial derivatives of ej (=normalized ebar_j) by using the derivatives of ebar_j
1154 for (int n=1; n<=dejdp.Length(); n++)
1155 {
1156 if (define_theta_by_a_transposed)
1157 {
1158 dejdp(n)(i) = deidp(n)(j);
1159 dejdp(n)(k) = dekdp(n)(j);
1160 }
1161
1162 //dejdp(n) = dejdp(n)*(1./ej_norm) - ej*((ej*dejdp(n))/ej_norm3);
1163 dejdp(n) *= 1./ej_norm;
1164 dejdp(n) -= ej*(ej*dejdp(n));
1165 }
1166
1167 if (!derivatives_wrt_rotated_vars)
1168 {
1169 //assembling of a help-vector saving at the (m*Dim()+n)th position the value
1170 //(ek^T * dejdpy(n))*SFy(m) + (ek^T * dejdpz(n))*SFz(m)
1171 for (int n=1; n<=Dim(); n++)
1172 {
1173 double fac_y = ek*dejdp(n);
1174 double fac_z = ek*dejdp(n+Dim());
1175 for (int m=1; m<=NS(); m++)
1176 {
1177 helpy((m-1)*Dim() + n) = fac_y*GetSFy(m, ploc) + fac_z*GetSFz(m, ploc);
1178 }
1179 }
1180 }
1181 else
1182 {
1183 //assembling of a help-vector saving at the (m*Dim()+n)th position the value
1184 //(ek^T * dejdpy(n))*A^T*SFy(m) + (ek^T * dejdpz(n))*A^T*SFz(m)
1185
1186 Vector3D vec_y, vec_z;
1187 for (int n=1; n<=Dim(); n++)
1188 {
1189 vec_y(n) = ek*dejdp(n);
1190 vec_z(n) = ek*dejdp(n+Dim());
1191 }
1192
1193 Matrix3D A = GetRotMatrix(ploc);
1194
1195 vec_y = vec_y*A.GetTp();
1196 vec_z = vec_z*A.GetTp();
1197
1198
1199 for (int n=1; n<=Dim(); n++)
1200 {
1201 for (int m=1; m<=NS(); m++)
1202 {
1203 helpy((m-1)*Dim() + n) = vec_y(n)*GetSFy(m, ploc) + vec_z(n)*GetSFz(m, ploc);
1204 }
1205 }
1206 }
1207
1208 //dimension of output matrix d is SOS() times Dim() (transpose of dAdq)
1209 for (int m=1; m<=SOS(); m++)
1210 {
1211 for (int n=1; n<=Dim(); n++)
1212 {
1213 d(m,n) += helpy(m)*ei(n);
1214 }
1215 }
1216 }
1217
1218
1219 ////-----------------------------------------
1220 ////old version:
1221 ////-----------------------------------------
1222 ////
1223 //Matrix3D A = GetRotMatrix(ploc);
1224 ////Matrix3D I = A.GetTp()*A;
1225 ////mbs->UO() << I.DoubleContract(I)-3. << "\n";
1226
1227 //Matrix drotvdqT(SOS(),Dim());
1228 //for (int i=1; i<=Dim(); i++)//i=1,2,3
1229 //{
1230 // int j = i%3 + 1;//i=1->j=2->k=3, i=2->j=3->k=1, i=3->j=1->k=2
1231 // int k = j%3 + 1;
1232 // Vector3D ei(A(1,i), A(2,i), A(3,i)); //i-th column of A
1233 // GetdRotvdqT(Vector3D(A(1,k), A(2,k), A(3,k)), ploc, drotvdqT, 1); //k-th column of A (=v) * dRotdq -> dRotvdq
1234 // for (int l=1; l<=SOS(); l++)//SOS=Dim*NS=3*9=27
1235 // for (int m=1; m<=Dim(); m++)//3
1236 // d(l,m) += drotvdqT(l,j)*ei(m);
1237 //}
1238 }
1239
1240
1241 //----------------
1242 //Massmatrix M
1243 //----------------
EvalM(Matrix & m,double t)1244 void ANCFBeamShear3DGeneric::EvalM(Matrix& m, double t)
1245 {
1246 if (massmatrix.Getcols() == SOS())
1247 {
1248 m = massmatrix;
1249 return;
1250 }
1251 else
1252 {
1253 int dim = Dim();
1254 int ns = NS();
1255 Vector SV;
1256 SV.SetLen(ns);
1257 Matrix HL(ns*dim,dim);
1258 Vector x1,x2,x3,w1,w2,w3; //integration points and according weights
1259
1260 //GetIntegrationRule => Gau� points on interval [-1,+1]
1261 GetIntegrationRule(x1,w1,7); //optimal: 6x3x3, 6x1x1 geht auch; (GetIntegrationRule-Order == 6 complies with == 7)
1262 GetIntegrationRule(x2,w2,3);
1263 GetIntegrationRule(x3,w3,3);
1264
1265 for (int i1=1; i1<=x1.GetLen(); i1++)
1266 {
1267 for (int i2=1; i2<=x2.GetLen(); i2++)
1268 {
1269 for (int i3=1; i3<=x3.GetLen(); i3++)
1270 {
1271 double x = x1(i1)*0.5*GetLx(); //x1 in [-1,+1], scaled rect. reference element: \xi
1272 double y = x2(i2)*0.5*GetLy();
1273 double z = x3(i3)*0.5*GetLz();
1274 Vector3D ploc(x,y,z);
1275 Vector3D rx0 = GetInitPosx3D(ploc);
1276 double rxn = rx0.Norm(); //curvature in reference element
1277 double jacdet = 0.5*GetLx()*0.5*GetLy()*0.5*GetLz() * rxn; //element transformation
1278
1279 GetShapes(SV,ploc);
1280 for (int i=0; i<ns; i++)
1281 {
1282 for (int j=1; j<=dim; j++)
1283 {
1284 HL(i*dim+j,j)=SV(i+1);
1285 }
1286 }
1287 m += fabs (jacdet) * Rho() * w1(i1)*w2(i2)*w3(i3) * (HL*HL.GetTp());
1288 }
1289 }
1290 }
1291 }
1292
1293 };
1294
GetJacobi(Matrix3D & jac,const Vector3D & ploc,const Vector & xg0) const1295 void ANCFBeamShear3DGeneric::GetJacobi(Matrix3D& jac, const Vector3D& ploc, const Vector& xg0) const
1296 {
1297 jac.SetSize(3,3);
1298 jac.SetAll(0.);
1299 for (int i = 1; i <= Dim(); i++)
1300 {
1301 for (int k=1; k <= NS(); k++)
1302 {
1303 jac(i,1) += GetSFx(k, ploc)*xg0((k-1)*Dim()+i);
1304 jac(i,2) += GetSFy(k, ploc)*xg0((k-1)*Dim()+i);
1305 jac(i,3) += GetSFz(k, ploc)*xg0((k-1)*Dim()+i);
1306 }
1307 }
1308 }
1309
EvalF2(Vector & f,double t)1310 void ANCFBeamShear3DGeneric::EvalF2(Vector& f, double t)
1311 {
1312 Body3D::EvalF2(f,t);
1313 //TMStartTimer(22); //CPU timing
1314 Vector f_elastic_forces(SOS()); //element residual - components are set to 0
1315
1316 if (InContinuumMechanicsFormulation())
1317 {
1318 EvalF2CM(f_elastic_forces, t);
1319 }
1320 else
1321 {
1322 EvalF2SM(f_elastic_forces, t);
1323 }
1324
1325 f -= f_elastic_forces; //f=f-fadd, (Ku=f -> Residual=f-Ku, Ku=fadd)
1326 //TMStopTimer(22);
1327 }
1328
1329 //residual of variational formulation in structural mechanics formulation
EvalF2SM(Vector & fadd,double t)1330 void ANCFBeamShear3DGeneric::EvalF2SM(Vector& fadd, double t)
1331 {
1332 if (!InContinuumMechanicsFormulation())
1333 //Reissner-Simo-VuQuoc Formulation
1334 {
1335
1336 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1337 //integration rule must be same as in GetFieldVariableValue(...) !!!! JG2013-09
1338 int order_x_bend = 3; //3
1339 int order_x_thickness_eps = 4; //4
1340 int order_x_axial_eps = 3; //3
1341 int order_x_shear = 3; //3
1342 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1343
1344 double beamGJ = GetMaterial().BeamGJkx();
1345 double beamEIy = GetMaterial().BeamEIy();
1346 double beamEIz = GetMaterial().BeamEIz();
1347 double beamEA = GetMaterial().BeamEA();
1348 double beamGAky = GetMaterial().BeamGAky();
1349 double beamGAkz = GetMaterial().BeamGAkz();
1350 // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1351 double beamGAkyz = BeamGAkyz(beamGAky, beamGAkz); // ACHTUNG: Literaturstudie notwendig - dieser Faktor wird f�r den letzten Term bzgl. Dickendeformation verwendet (GAkyz*2E_yz*2deltaE_yz).
1352 // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1353
1354 //UO() << "beamEIy=" << beamEIy << "\n";
1355 //UO() << "beamEIz=" << beamEIz << "\n";
1356
1357 Vector faddadd(SOS());
1358 Vector3D gk(0.); //container for Gamma, Kappa, and Strain E=(E_yy,E_zz,E_yz) respectively
1359 Matrix delta_gk(SOS(),Dim()); //container for DeltaGamma, and DeltaKappa, and DeltaStrain respectively
1360 Vector wT,xT; //integration points and weights
1361 GetIntegrationRule(xT,wT,order_x_bend);
1362
1363 // add kappa terms (torsion + 2*bending)
1364 //TMStartTimer(23); //CPU timing %52%
1365 for (int i1=1; i1<=xT.GetLen(); i1++) //i1-loop over all integration points
1366 {
1367 double x = xT(i1); //x in unit configuration [-1 , 1]
1368 double xi = x*GetLx()*0.5; //xi in scaled configuration [-lx/2 , lx/2]
1369 double det = 0.5*GetLx(); //determinant of the jacobian (from element transformation [-1 , 1] -> [-lx/2 , lx/2])
1370
1371 double fact1 = penalty(1)*beamGJ*wT(i1)*det; // beamGJ = torsional stiffness (including torsional correction factor, depending on the shape of the cross section)
1372 double fact2 = penalty(2)*beamEIy*wT(i1)*det; // beamEIy = bending stiffness (in e2=y direction))
1373 double fact3 = penalty(3)*beamEIz*wT(i1)*det; // beamEIz = bending stiffness (in e3=z direction)
1374
1375 GetDeltaKappa(xi, delta_gk, gk);
1376 //mbs->UO() << "int-point" << i1 << "\n";
1377 //if(i1 == 2)
1378 //{
1379 // //mbs->UO() << "test\n";
1380 // ofstream ofs; ofs.open("..\\..\\output\\FWF\\kappa_bei_GC12_Biegung_um_y_linelem_ip3.txt");
1381 // ofs << "gk = " << gk << "\n" << "delta_gk = " << delta_gk << "\n";
1382 // ofs.close();
1383 //}
1384 gk.Scale(1./fact1, 1./fact2, 1./fact3);
1385
1386 Mult(delta_gk, gk, faddadd);
1387 fadd += faddadd;
1388 }
1389
1390 // add gamma terms (axial strains only)
1391 GetIntegrationRule(xT,wT,order_x_axial_eps);
1392 //GetIntegrationRuleLobatto(xT,wT,2);
1393
1394 for (int i1=1; i1<=xT.GetLen(); i1++) //i1-loop over all integration points
1395 {
1396 double x = xT(i1); //x in unit configuration [-1 , 1]
1397 double xi = x*GetLx()*0.5; //xi in scaled configuration [-lx/2 , lx/2]
1398 double det = 0.5*GetLx(); //determinant of the jacobian (from element transformation [-1 , 1] -> [-lx/2 , lx/2])
1399
1400 double fact1 = penalty(4)*beamEA*wT(i1)*det; // beamEA = axial stiffness
1401 //integration of axial term only!!!
1402 double fact2 = 0.; //Gamma2 and Gamma3 is integrated seperately
1403 double fact3 = 0.;
1404
1405 //GetGamma(xi, gk);
1406 GetDeltaGamma(xi, delta_gk, gk);
1407 gk.Scale(1./fact1, 1./fact2, 1./fact3);
1408 //UO() << "gk=" << gk << "\n";
1409 Mult(delta_gk, gk, faddadd);
1410
1411 fadd += faddadd;
1412 }
1413
1414 // add gamma terms (shear terms only)
1415 GetIntegrationRule(xT,wT,order_x_shear);
1416 //GetIntegrationRuleLobatto(xT,wT,order_x_shear);
1417 for (int i1=1; i1<=xT.GetLen(); i1++) //i1-loop over all integration points
1418 {
1419 double x = xT(i1); //x in unit configuration [-1 , 1]
1420 double xi = x*GetLx()*0.5; //xi in scaled configuration [-lx/2 , lx/2]
1421 double det = 0.5*GetLx(); //determinant of the jacobian (from element transformation [-1 , 1] -> [-lx/2 , lx/2])
1422
1423 double fact1 = 0.; //integration of shear terms only!!!
1424 double fact2 = penalty(5)*beamGAky*wT(i1)*det;
1425 double fact3 = penalty(6)*beamGAkz*wT(i1)*det;
1426
1427 GetDeltaGamma(xi, delta_gk, gk);
1428 gk.Scale(1./fact1, 1./fact2, 1./fact3);
1429 //UO() << "gk=" << gk << "\n";
1430 Mult(delta_gk, gk, faddadd);
1431 fadd += faddadd;
1432 }
1433 //TMStopTimer(24); //CPU timing
1434
1435 // add thickness terms (E_yy E_zz 2*E_yz)
1436 //TMStartTimer(25); //CPU timing %7.5%
1437 //GetIntegrationRule(xT,wT,order_x_thickness_eps);
1438 GetIntegrationRuleLobatto(xT,wT,order_x_thickness_eps);
1439 for (int i1=1; i1<=xT.GetLen(); i1++) //i1-loop over all integration points
1440 {
1441 double x = xT(i1); //x in unit configuration [-1 , 1]
1442 double xi = x*GetLx()*0.5; //xi in scaled configuration [-lx/2 , lx/2]
1443 double det = 0.5*GetLx(); //determinant of the jacobian (from element transformation [-1 , 1] -> [-lx/2 , lx/2])
1444
1445 double fact1 = penalty(7)*beamEA*wT(i1)*det;
1446 double fact2 = penalty(8)*beamEA*wT(i1)*det;
1447 double fact3 = penalty(9)*beamGAkyz*wT(i1)*det; // ACHTUNG: beamGAkyz: Literaturstudie noch notwendig!!!!!!!
1448
1449 GetDeltaThicknessStrain(xi, delta_gk, gk);
1450 gk.Scale(1./fact1, 1./fact2, 1./fact3);
1451
1452 Mult(delta_gk, gk, faddadd);
1453 fadd += faddadd;
1454 }
1455 //TMStopTimer(25); //CPU timing
1456 }
1457 }
1458
1459 //residual of variational formulation in continuum mechanics formulation
EvalF2CM(Vector & fadd,double t)1460 void ANCFBeamShear3DGeneric::EvalF2CM(Vector& fadd, double t)
1461 {
1462 int dim = Dim();//=3
1463 int ns = NS();//=9
1464 int sos = SOS();//9*3
1465
1466 ConstVector<27> temp;//KN-todo: 27 ersetzen durch ConstInt MaxDOF ...
1467 temp.SetLen(SOS());
1468 temp.SetAll(0);
1469
1470 double A = GetMaterial().BeamA(); //cross section area
1471 double EA = GetMaterial().BeamEA(); //Em*area
1472 double GAky = GetMaterial().BeamGAky();
1473 double GAkz = GetMaterial().BeamGAkz();
1474 double G = Em()/(2.+2.*Nu());
1475
1476 SetComputeCoordinates();//xg(i) = XG(i);
1477
1478 Matrix3D strain, piola1, F;
1479 strain.SetAll(0);
1480 piola1.SetAll(0);
1481
1482 int poissoncorrection = 0; //1==reduced integration of poisson part
1483 if (BeamFormulation == 2) poissoncorrection = 1;
1484
1485 static Vector x1,x2,x3,w1,w2,w3;//KN-todo: ersetzen durch ConstVector
1486
1487 //Vector3D ks(1.,GAky/(G*A),GAkz/(G*A)); //vector of shear correction factors
1488 //Matrix3D R = GetRotMatrix(Vector3D(0.,0.,0.), 1).GetTp(); //WARNING:for curved beams, this must be evaluated at every point of integration!
1489 //Vector3D ksR = R*ks;
1490 //ksR.X() = fabs(ksR.X());
1491 //ksR.Y() = fabs(ksR.Y());
1492 //ksR.Z() = fabs(ksR.Z());
1493 //UO() << "rot=" << R << "\n";
1494 //UO() << "ks=" << ks << "\n";
1495 //UO() << "ksR=" << ksR << "\n";
1496
1497 for (int kk=1; kk <= 1+poissoncorrection; kk++)
1498 {
1499 ConstMatrix<36> Dm;
1500 Dm.SetSize(6,6);
1501 Dm.SetAll(0.);
1502 if (poissoncorrection)
1503 {
1504 if (BeamFormulation == 2)//locking compensated
1505 {
1506 if (kk == 1)//D^0 (equ. (25))
1507 {
1508 //poisson ratio zero:
1509 Dm(1,1) = Em();
1510 Dm(2,2) = Em();//reduced thickness stiffness 0.5*
1511 Dm(3,3) = Em();//reduced thickness stiffness
1512
1513 //Dm(4,4) = G*5./6.;//1.;//*ksR.X();//statt GAks/A; //gamma_yz
1514 //Dm(5,5) = G*5./6.;//*ksR.Y(); //gamma_xz
1515 //Dm(6,6) = G*5./6.;//*ksR.Z();//event. weglassen -> konstant �ber Querschnitt //gamma_xy
1516 //UO() << "CMF (locking-free)" << "\n";
1517 Dm(4,4) = G;//statt GAks/A; //gamma_yz
1518 Dm(5,5) = GAky/A; //gamma_xz
1519 Dm(6,6) = GAkz/A;//event. weglassen -> konstant �ber Querschnitt //gamma_xy
1520
1521 //x_i ... integration points, w_i ... integration weights
1522 GetIntegrationRule(x1,w1,order_axial);//for dx
1523 GetIntegrationRule(x2,w2,order_crosssectional);//for dy
1524 GetIntegrationRule(x3,w3,order_crosssectional);//for dz
1525 }
1526 else if (kk == 2)//D^v (equ. (26))
1527 {
1528 //$ KN 2011-10-28: help1, help2 mit Paper verglichen -> passt
1529 double help1 = (Em()* (1-Nu()) / ((1+Nu())*(1-2.*Nu()))) - Em();
1530 //help1 = (2.*Em()*Nu()*Nu())/((1+Nu())*(1-2.*Nu()))) ... equ.(26)
1531 Dm(1,1) = help1;
1532 Dm(2,2) = help1;
1533 Dm(3,3) = help1;
1534
1535 double help2 = Em()*Nu()/((1+Nu())*(1-2.*Nu()));
1536 Dm(1,2) = help2;
1537 Dm(2,1) = help2;
1538 Dm(1,3) = help2;
1539 Dm(3,1) = help2;
1540 Dm(2,3) = help2;
1541 Dm(3,2) = help2;
1542
1543 GetIntegrationRule(x1,w1,order_axial);//x_i ... integration points, w_i ... integration weights
1544 GetIntegrationRule(x2,w2,1);//only 1 int. point along cross section//Gau�Integration
1545 GetIntegrationRule(x3,w3,1);
1546 }
1547 }
1548 }
1549 else if(BeamFormulation == 1)//cont. mech. approach which shows locking phenomenon
1550 {
1551 //x_i ... integration points, w_i ... integration weights
1552 GetIntegrationRule(x1,w1,order_axial);//for dx
1553 GetIntegrationRule(x2,w2,order_crosssectional);//for dy//Gau�!!!
1554 GetIntegrationRule(x3,w3,order_crosssectional);//for dz
1555
1556 double help3 = Em()* (1-Nu()) / ((1+Nu())*(1-2*Nu()));
1557 double help4 = Em()*Nu() / ((1+Nu())*(1-2*Nu()));
1558
1559 Dm(1,1)=help3;
1560 Dm(1,2)=help4;
1561 Dm(1,3)=help4;
1562
1563 Dm(2,1)=help4;
1564 Dm(2,2)=help3;
1565 Dm(2,3)=help4;
1566
1567 Dm(3,1)=help4;
1568 Dm(3,2)=help4;
1569 Dm(3,3)=help3;
1570
1571 Dm(4,4) = G;//statt GAks/A; //gamma_yz
1572 Dm(5,5) = GAky/A; //gamma_xz
1573 Dm(6,6) = GAkz/A;//event. weglassen -> konstant �ber Querschnitt //gamma_xy
1574 /*UO() << "false!\n";*/
1575 //Dm(4,4)=G*5./6.;
1576 //Dm(5,5)=G*5./6.;
1577 //Dm(6,6)=G*5./6.;
1578 //Dm(4,4)=10*Em();//G*5./6.;
1579 //Dm(5,5)=10*Em();//G*5./6.;
1580 //Dm(6,6)=10*Em();//G*5./6.;
1581 //UO() << "D=" << Dm << "\n";
1582 }
1583 else if(BeamFormulation == 3)// only D^0 (equ. (25))
1584 {
1585 Dm(1,1) = Em();
1586 Dm(2,2) = Em();
1587 Dm(3,3) = Em();
1588 //Dm(4,4) = G*ksR.X();//statt GAks/A; //gamma_yz
1589 //Dm(5,5) = G*ksR.Y(); //gamma_xz
1590 //Dm(6,6) = G*ksR.Z();//event. weglassen -> konstant �ber Querschnitt //gamma_xy
1591 /*UO() << "false!\n";*/
1592 Dm(4,4) = G*5./6.;//statt GAks/A;
1593 Dm(5,5) = G*5./6.;
1594 Dm(6,6) = G*5./6.;
1595
1596 //x_i ... integration points, w_i ... integration weights
1597 GetIntegrationRule(x1,w1,order_axial);//for dx
1598 GetIntegrationRule(x2,w2,order_crosssectional);//for dy
1599 GetIntegrationRule(x3,w3,order_crosssectional);//for dz
1600 }
1601 //mbs->UO()<<"Dm="<<Dm<<"\n";
1602
1603 for (int i1=1; i1<=x1.GetLen(); i1++)
1604 {
1605 for (int i2=1; i2<=x2.GetLen(); i2++)
1606 {
1607 for ( int i3=1; i3<=x3.GetLen(); i3++)
1608 {
1609 int i,j,k;
1610 Vector3D p(x1(i1)*0.5*GetLx(), x2(i2)*0.5*GetLy(), x3(i3)*0.5*GetLz());
1611
1612 Matrix3D jac;
1613 GetJacobi(jac,p,q0);//p(1) in [-lx/2,lx/2]
1614 double jacdet = jac.Det();
1615 //mbs->UO()<<"jacdet = "<<jacdet<<"\n";
1616 //mbs->UO()<<"0.5*lx*0.5*ly*0.5*lz = "<<0.5*lx*0.5*ly*0.5*lz<<"\n";
1617
1618 Matrix3D jacinv;
1619 jac.GetInverse(jacinv);
1620 jacinv = jacinv.GetTp();
1621 //mbs->UO() << "jacinv = " << jacinv << "\n";
1622
1623 static Matrix grad0, grad;//am Einheitselement
1624 grad0.SetSize(Dim(),NS());
1625 for(i=1; i<=NS(); i++)
1626 {
1627 grad0(1,i)=GetSFx(i,p);//3x9-matrix
1628 grad0(2,i)=GetSFy(i,p);
1629 grad0(3,i)=GetSFz(i,p);
1630 }
1631
1632 Mult(jacinv, grad0, grad);//grad0*jacinv
1633 F.SetAll(0);//muss in jedem Integrationsschritt zur�ck auf 0 gesetzt werden!!!!
1634 for (int j = 1; j <= dim; j++)
1635 {
1636 for (int i = 1; i <= NS(); i++)
1637 {
1638 int l = (i-1)*dim+j;
1639 F(j,1) += grad(1,i)*xg(l);//3x9-matrix*9-vector
1640 F(j,2) += grad(2,i)*xg(l);
1641 F(j,3) += grad(3,i)*xg(l);
1642 }
1643 }
1644 //F = F.GetTp();//FTranspose?
1645 if(calclinear == 0)
1646 {
1647 F(1,1) += 1;//Positionsgrad (nicht Grad(u), sondern Grad(r))
1648 F(2,2) += 1;
1649 F(3,3) += 1;
1650 //E = 1/2*(F^T F-I)
1651 strain = 0.5*(F.GetTp()*F);//eps=0.5*(F.GetTp()*F-Id)
1652 strain(1,1) -= 0.5;
1653 strain(2,2) -= 0.5;
1654 strain(3,3) -= 0.5;
1655 }
1656 else if(calclinear == 1)
1657 {
1658 strain = 0.5*(F + F.GetTp());
1659 }
1660
1661 ConstVector<6> s3(strain(1,1), strain(2,2), strain(3,3), 2.*strain(2,3), 2.*strain(1,3), 2.*strain(1,2));
1662 ConstVector<6> stress;
1663 stress = Dm * s3;//Dm=6x6-Matrix, //S = D : E
1664
1665 //generate 2. Piola-Kirchhoff-tensor, 2.PK = symm. tensor
1666 piola1(1,1)=stress(1);
1667 piola1(2,2)=stress(2);
1668 piola1(3,3)=stress(3);
1669
1670 piola1(2,3)=stress(4);
1671 piola1(3,2)=stress(4);
1672 piola1(1,3)=stress(5);
1673 piola1(3,1)=stress(5);
1674 piola1(1,2)=stress(6);
1675 piola1(2,1)=stress(6);
1676
1677 if(calclinear == 0)
1678 {
1679 piola1 = F * piola1;//P = F * S, 1.PK = F * 2.PK
1680 }
1681
1682 //grad = grad.GetTp();//gradTranspose?
1683 for (int j=1; j <= dim; j++)
1684 {
1685 for (int i = 0; i < ns; i++)
1686 {
1687 temp(dim*i+j) = grad(1,i+1)*piola1(j,1) + grad(2,i+1)*piola1(j,2) + grad(3,i+1)*piola1(j,3);
1688 }
1689 }
1690 //integrate(P : dF/dq)
1691 fadd.MultAdd(fabs(jacdet) * 0.5*GetLx()*0.5*GetLy()*0.5*GetLz() * w1(i1)*w2(i2)*w3(i3),temp);
1692
1693 }
1694 }
1695 }
1696 }
1697 };
1698
1699 //----------------------------------------------------------------------------
1700 //Gamma & Kappa & ThicknessStrain + Derivatives of those with respect to dof
1701 //----------------------------------------------------------------------------
GetGamma(const double & xi,Vector3D & gamma,int flagD)1702 void ANCFBeamShear3DGeneric::GetGamma(const double& xi, Vector3D& gamma, int flagD) // -Lx/2 <= xi <= Lx/2
1703 {
1704 //gamma1 = t1*drdx/|e1| - 1
1705 //gamma2 = t2*drdxi
1706 //gamma3 = t3*drdxi
1707 //
1708 //gamma = (gamma1, gamma2, gamma3)
1709 Vector3D p0(xi, 0., 0.); //xi in [-Lx/2, Lx/2]
1710
1711 Vector3D drdx, drdy, drdz, e1, e2, e3;
1712 drdx = GetPosx3D(p0, flagD);
1713 drdy = GetPosy3D(p0, flagD);
1714 drdz = GetPosz3D(p0, flagD);
1715 e1 = Getei(1,
1716 drdy.X(), drdy.Y(), drdy.Z(),
1717 drdz.X(), drdz.Y(), drdz.Z());
1718 e2 = Getei(2,
1719 drdy.X(), drdy.Y(), drdy.Z(),
1720 drdz.X(), drdz.Y(), drdz.Z());
1721 e3 = Getei(3,
1722 drdy.X(), drdy.Y(), drdy.Z(),
1723 drdz.X(), drdz.Y(), drdz.Z());
1724 e1.Normalize();
1725 e2.Normalize();
1726 e3.Normalize();
1727
1728 gamma.Set(
1729 e1*drdx - 1.,
1730 e2*drdx,
1731 e3*drdx);
1732 }
1733
GetDeltaGamma(const double & xi,Matrix & dgammadq,Vector3D & gamma)1734 void ANCFBeamShear3DGeneric::GetDeltaGamma(const double& xi, Matrix& dgammadq, Vector3D& gamma) // -Lx/2 <= xi <= Lx/2
1735 {
1736 //dgamma1dq = ( ... ) see paper
1737 //dgamma2dq = ( r,y.S,x + r,x.S,y - (r,x.r,y)/(r,y.r,y).r,y.S,y ) / (norm(r,y))
1738 //dgamma3dq = ( r,z.S,x + r,x.S,z - (r,x.r,z)/(r,z.r,z).r,z.S,z ) / (norm(r,z))
1739 //
1740 //dgammadq = (dgamma1dq, dgamma2dq, dgamma3dq)
1741
1742 Vector3D p0(xi, 0., 0.); //xi in [-Lx/2, Lx/2]
1743 Vector sfx(NS(), 1); GetShapesx(sfx, p0);
1744 Vector sfy(NS(), 1); GetShapesy(sfy, p0);
1745 Vector sfz(NS(), 1); GetShapesz(sfz, p0);
1746
1747 dgammadq.SetSize(SOS(), Dim());
1748 dgammadq.SetAll(0.);
1749
1750 Vector3D drdx, drdy, drdz;
1751 drdx = GetPosx3D(p0);
1752 drdy = GetPosy3D(p0);
1753 drdz = GetPosz3D(p0);
1754
1755 // help variables
1756 double norm_drdy_squ = drdy.Norm2();
1757 double norm_drdy_lin = sqrt(norm_drdy_squ);
1758 double norm_drdy_cub = norm_drdy_lin * norm_drdy_squ;
1759 double norm_drdz_squ = drdz.Norm2();
1760 double norm_drdz_lin = sqrt(norm_drdz_squ);
1761 double norm_drdz_cub = norm_drdz_lin * norm_drdz_squ;
1762
1763 for (int i=1; i<=Dim(); i++)
1764 {
1765 Vector3D ei = Getei(i, drdy.X(), drdy.Y(), drdy.Z(), drdz.X(), drdz.Y(), drdz.Z());
1766 TArray<Vector3D> deidp(GetdeidDiffVar(i, drdy.X(), drdy.Y(), drdy.Z(), drdz.X(), drdz.Y(), drdz.Z()));
1767 double norm_ei_squ = ei.Norm2();
1768 double norm_ei_lin = sqrt(norm_ei_squ);
1769 double norm_ei_cub = norm_ei_lin * norm_ei_squ;
1770
1771 gamma(i) = (ei*drdx)/norm_ei_lin;
1772 Vector3D v = (1/norm_ei_lin)*drdx - (gamma(i)/norm_ei_squ)*ei;
1773 //double facy = v * (de1dp(1) + de1dp(2) + de1dp(3));
1774 //double facz = v * (de1dp(4) + de1dp(5) + de1dp(6));
1775
1776 //assembling of dgamma1dq
1777 for (int j=1; j<=NS(); j++)
1778 {
1779 for (int k=1; k<=Dim(); k++)
1780 {
1781 double rowidx = (j-1)*Dim() + k;
1782 dgammadq(rowidx, i) = ((1/norm_ei_lin)*ei(k))*sfx(j);
1783 dgammadq(rowidx, i) += (v*deidp(k))*sfy(j);
1784 dgammadq(rowidx, i) += (v*deidp(k+Dim()))*sfz(j);
1785 }
1786 }
1787 }
1788 gamma(1) -= 1;
1789 }
1790
1791
GetKappa(const double & xi,Vector3D & kappa,int flagD)1792 void ANCFBeamShear3DGeneric::GetKappa(const double& xi, Vector3D& kappa, int flagD) // -Lx/2 <= xi <= Lx/2
1793 {
1794 Vector3D p0 = Vector3D(xi, 0., 0.);
1795 kappa = GetTwistAndCurvature(p0, 0, flagD) - GetTwistAndCurvature(p0, 1, flagD);
1796 }
1797
1798
GetDeltaKappa(const double & xi,Matrix & dkappadq,Vector3D & kappa)1799 void ANCFBeamShear3DGeneric::GetDeltaKappa(const double& xi, Matrix& dkappadq, Vector3D& kappa) //derivative of vector Kappa wrt vector q at xi (which is between -Lx/2 and Lx/2)
1800 {
1801 //TMStartTimer(28);
1802 Vector3D p0 = Vector3D(xi, 0, 0);
1803
1804 Vector sfy(NS(), 1); GetShapesy(sfy, p0);
1805 Vector sfz(NS(), 1); GetShapesz(sfz, p0);
1806 Vector sfxy(NS(), 1); GetShapesxy(sfxy, p0);
1807 Vector sfxz(NS(), 1); GetShapesxz(sfxz, p0);
1808
1809
1810 Vector3D drdy = GetPosy3D(p0);
1811 Vector3D drdz = GetPosz3D(p0);
1812 Vector3D ddrdxdy = GetPosxy3D(p0);
1813 Vector3D ddrdxdz = GetPosxz3D(p0);
1814
1815 // e1 = drdy x drdz, e3 = drdz, e2 = e3 x e1 = -e1 x e3
1816 Vector3D ei;
1817 Vector3D deidxi;
1818 TArray<Vector3D> deidp(6);
1819 TArray<Vector3D> ddeidxidp(12);
1820
1821 // ei depends on p = (dr1dy dr2dy dr3dy dr1dz dr2dz dr3dz),
1822 // deidxi depends on p = (dr1dy dr2dy dr3dy dr1dz dr2dz dr3dz ddr1dxdy ddr2dxdy ddr3dxdy ddr1dxdz ddr2dxdz ddr3dxdz)
1823 // delta kappa = delta k - delta(Rot*k_0)
1824
1825 // calculation of delta k: in the way dkdq = dkdp * dpdq
1826
1827 TArray<Vector3D> dkdp(12); // Vector3D() sets components to 0!
1828 for (int i=1; i<=Dim(); i++)
1829 {
1830 ei = Getei(i,
1831 drdy.X(), drdy.Y(), drdy.Z(),
1832 drdz.X(), drdz.Y(), drdz.Z());
1833 deidxi = Getdeidxi(i,
1834 drdy.X(), drdy.Y(), drdy.Z(),
1835 drdz.X(), drdz.Y(), drdz.Z(),
1836 ddrdxdy.X(), ddrdxdy.Y(), ddrdxdy.Z(),
1837 ddrdxdz.X(), ddrdxdz.Y(), ddrdxdz.Z());
1838 deidp = GetdeidDiffVar(i,
1839 drdy.X(), drdy.Y(), drdy.Z(),
1840 drdz.X(), drdz.Y(), drdz.Z());
1841 ddeidxidp = GetddeidxidDiffVar(i,
1842 drdy.X(), drdy.Y(), drdy.Z(),
1843 drdz.X(), drdz.Y(), drdz.Z(),
1844 ddrdxdy.X(), ddrdxdy.Y(), ddrdxdy.Z(),
1845 ddrdxdz.X(), ddrdxdz.Y(), ddrdxdz.Z());
1846
1847 double einormsquare = ei*ei;
1848
1849 for (int j=1; j<=12; j++)
1850 {
1851 if (j<=6)
1852 {
1853 dkdp(j) -= ((ei*deidp(j)) / Sqr(einormsquare)) * ei.Cross(deidxi);
1854 dkdp(j) += (0.5/einormsquare) * deidp(j).Cross(deidxi);
1855 }
1856
1857 dkdp(j) += (0.5/einormsquare) * ei.Cross(ddeidxidp(j));
1858 }
1859 }
1860
1861 // delta k = [ d k / d p ] . [ d p / dq ]
1862 dkappadq.SetAll(0.);
1863 Vector3D rowvec;
1864
1865 for(int k=1; k<=NS(); k++)
1866 {
1867 for(int j=1; j<=Dim(); j++)
1868 {
1869 rowvec =
1870 dkdp(j )*sfy(k) +
1871 dkdp(j+3)*sfz(k) +
1872 dkdp(j+6)*sfxy(k) +
1873 dkdp(j+9)*sfxz(k);
1874 int rowidx = j + (k-1)*Dim();
1875
1876 dkappadq(rowidx, 1) = rowvec.X();
1877 dkappadq(rowidx, 2) = rowvec.Y();
1878 dkappadq(rowidx, 3) = rowvec.Z();
1879 }
1880 }
1881
1882 // [kappa]_e = A^T (k - k_0) [..] be the representation in the local frame basis
1883 // d [kappa]_e / d q = d A^T /dq (k - k_0) + A^T dk/dq
1884
1885 Matrix helper = dkappadq;
1886 Matrix3D A = GetRotMatrix(p0);
1887 Mult(helper, A, dkappadq); //dkappadq = helper * A;
1888
1889 //Matrix3D AT = A;//.GetTp();
1890 //for(int i=1; i<=NS()*Dim(); i++)
1891 //{
1892 // Vector3D v;
1893 // v(1)=helper(i,1);
1894 // v(2)=helper(i,2);
1895 // v(3)=helper(i,3);
1896 // v = AT*v;
1897 // dkappadq(i,1) = v(1);
1898 // dkappadq(i,2) = v(2);
1899 // dkappadq(i,3) = v(3);
1900 //}
1901
1902 //TMStopTimer(29);
1903
1904 //TMStartTimer(30);
1905 kappa = GetTwistAndCurvature(p0, 0) - GetTwistAndCurvature(p0, 1);
1906
1907 //TMStopTimer(30);
1908 //TMStartTimer(20);
1909 GetdRotvdqT(A*kappa, p0, helper, 1);
1910 //TMStopTimer(20);
1911 dkappadq += helper;
1912
1913 //kappa = A.GetTp()*kappa;
1914 }
1915
1916
GetThicknessStrain(const double & xi,Vector3D & thicknessstrain,int flagD)1917 void ANCFBeamShear3DGeneric::GetThicknessStrain(const double& xi, Vector3D& thicknessstrain, int flagD) // -Lx/2 <= xi <= Lx/2
1918 {
1919 //thicknessstrain_yy = 1/2 * (drdy*drdy - 1)
1920 //thicknessstrain_zz = 1/2 * (drdz*drdz - 1)
1921 //thicknessstrain_yz_twice = 2 * thicknessstrain_yz = drdy*drdz
1922 //
1923 //thicknessstrain = (thicknessstrain_yy, thicknessstrain_zz, thicknessstrain_yz_twice)
1924
1925 Vector3D p0(xi, 0., 0.); //xi in [-Lx/2, Lx/2]
1926
1927 Vector3D drdy = GetPosy3D(p0, flagD);
1928 Vector3D drdz = GetPosz3D(p0, flagD);
1929
1930 thicknessstrain.Set(
1931 (drdy*drdy - 1.)/2.,
1932 (drdz*drdz - 1.)/2.,
1933 drdy*drdz );
1934 }
1935
GetDeltaThicknessStrain(const double & xi,Matrix & dthicknessstraindq,Vector3D & thicknessstrain)1936 void ANCFBeamShear3DGeneric::GetDeltaThicknessStrain(const double& xi, Matrix& dthicknessstraindq, Vector3D& thicknessstrain) // -Lx/2 <= xi <= Lx/2
1937 {
1938 //dthicknessstrainyydq = r,y.S,y
1939 //dthicknessstrainzzdq = r,z.S,z
1940 //dthicknessstrainyztwicedq = r,y.S,z + r,z.S,y
1941 //
1942 //dthicknessstraindq = (dthicknessstrainyydq, dthicknessstrainzzdq, dthicknessstrainyztwicedq)
1943
1944 Vector3D p0(xi, 0., 0.); //xi in [-Lx/2, Lx/2]
1945 Vector sfy(NS(), 1); GetShapesy(sfy, p0);
1946 Vector sfz(NS(), 1); GetShapesz(sfz, p0);
1947
1948 dthicknessstraindq.SetSize(SOS(), Dim());
1949 dthicknessstraindq.SetAll(0.);
1950
1951 Vector3D drdy = GetPosy3D(p0);
1952 Vector3D drdz = GetPosz3D(p0);
1953
1954 //set thickness strain
1955 thicknessstrain.Set(
1956 (drdy*drdy - 1.)/2.,
1957 (drdz*drdz - 1.)/2.,
1958 drdy*drdz );
1959
1960 //assembling of dthicknessstraindq
1961 for (int i=1; i<=NS(); i++)
1962 {
1963 for (int j=1; j<=Dim(); j++)
1964 {
1965 double rowidx = (i-1)*Dim() + j;
1966 dthicknessstraindq(rowidx, 1) += drdy(j)*sfy(i);
1967 dthicknessstraindq(rowidx, 2) += drdz(j)*sfz(i);
1968 dthicknessstraindq(rowidx, 3) += drdy(j)*sfz(i) + drdz(j)*sfy(i);
1969 }
1970 }
1971 }
1972
1973
1974 //for plasticity:
PostNewtonStep(double t)1975 double ANCFBeamShear3DGeneric::PostNewtonStep(double t) // changes plastic strains, returns yieldfunction(sigma)/Em
1976 {
1977 return 0;
1978 }
1979
1980
1981 //---------------------------------------------------------------
1982 //for Displacement/Stress/Strain/Stress Resultants-Visualization:
1983 //---------------------------------------------------------------
GetFieldVariableValue(const FieldVariableDescriptor & fvd,const Vector3D & local_position,bool flagD)1984 double ANCFBeamShear3DGeneric::GetFieldVariableValue(const FieldVariableDescriptor & fvd, const Vector3D & local_position, bool flagD)
1985 {
1986
1987 //ploc_scaled = local position: -l/2 <= ploc_scaled <= +l/2
1988 Vector3D p0 = local_position;//-l/2 <= p0 <= +l/2
1989
1990 xgd.SetLen(SOS());
1991 GetDrawCoordinates(xgd); //xgd=XGD
1992
1993 //Position:
1994 if (fvd.VariableType() == FieldVariableDescriptor::FVT_position)
1995 {
1996 Vector3D u;
1997 if(flagD)
1998 u = GetPosD(p0);
1999 else
2000 u = GetPos(p0);
2001
2002 return fvd.GetComponent(u);
2003 }
2004
2005 //Displacement:
2006 if (fvd.VariableType() == FieldVariableDescriptor::FVT_displacement)
2007 {
2008 Vector3D u;
2009 if(flagD)
2010 u = GetDisplacementD(p0);
2011 else
2012 u = GetDisplacement(p0);
2013
2014 return fvd.GetComponent(u);
2015 }
2016
2017 //Velocity:
2018 if (fvd.VariableType() == FieldVariableDescriptor::FVT_velocity)
2019 {
2020 Vector3D u;
2021 u = GetVel3D(p0, flagD);
2022
2023 return fvd.GetComponent(u);
2024 }
2025
2026 //Stress and Strain:
2027 if (InContinuumMechanicsFormulation())
2028 {
2029 Matrix3D strain(0), F, piola1;
2030 ConstMatrix<36> Dm;
2031 Dm.SetSize(6,6);
2032 Dm.SetAll(0);
2033 Vector3D rx, ry, rz;
2034 rx=GetPosx3D(p0,flagD);
2035 ry=GetPosy3D(p0,flagD);
2036 rz=GetPosz3D(p0,flagD);
2037
2038 //strain:
2039 F.SetAll(0); //$JG: WRONG: this is the deformation gradient in local coordinates
2040 for (int j = 1; j <= Dim(); j++)
2041 {
2042 F(j,1) = rx(j);//grad*xg
2043 F(j,2) = ry(j);
2044 F(j,3) = rz(j);
2045 }
2046 if(calclinear == 0)
2047 {
2048 //F(1,1) += 1;//Positionsgrad (nicht Grad(u), sondern Grad(r))//ux->rx: +1, hier haben wir schon rx, kein +1
2049 //F(2,2) += 1;
2050 //F(3,3) += 1;
2051 //E = 1/2*(F^T F-I)
2052 strain = 0.5*(F.GetTp()*F);//eps=0.5*(F.GetTp()*F-Id)
2053 strain(1,1) -= 0.5;
2054 strain(2,2) -= 0.5;
2055 strain(3,3) -= 0.5;
2056 }
2057 else if(calclinear == 1)
2058 {
2059 F(1,1) -= 1;//rx->ux: -1
2060 F(2,2) -= 1;
2061 F(3,3) -= 1;
2062 strain = 0.5*(F + F.GetTp());
2063 }
2064 //stress:
2065 double help3 = Em()* (1-Nu()) / ((1+Nu())*(1-2*Nu()));
2066 double help4 = Em()*Nu() / ((1+Nu())*(1-2*Nu()));
2067 double GAky = GetMaterial().BeamGAky();
2068 double GAkz = GetMaterial().BeamGAkz();
2069 double G = Em()/(2.+2.*Nu());
2070 double A = GetMaterial().BeamA();
2071
2072 Dm(1,1)=help3;
2073 Dm(1,2)=help4;
2074 Dm(1,3)=help4;
2075
2076 Dm(2,1)=help4;
2077 Dm(2,2)=help3;
2078 Dm(2,3)=help4;
2079
2080 Dm(3,1)=help4;
2081 Dm(3,2)=help4;
2082 Dm(3,3)=help3;
2083
2084 Dm(4,4)=G;
2085 Dm(5,5)=G;
2086 Dm(6,6)=G;
2087
2088 ConstVector<6> s3(strain(1,1), strain(2,2), strain(3,3), 2.*strain(2,3), 2.*strain(1,3), 2.*strain(1,2));
2089 ConstVector<6> stress;
2090 stress = Dm * s3;//Dm=6x6-Matrix, //S = D : E
2091
2092 //generate 2. Piola-Kirchhoff-tensor, 2.PK = symm. tensor
2093 piola1(1,1)=stress(1);
2094 piola1(2,2)=stress(2);
2095 piola1(3,3)=stress(3);
2096
2097 piola1(2,3)=stress(4);
2098 piola1(3,2)=stress(4);
2099 piola1(1,3)=stress(5);
2100 piola1(3,1)=stress(5);
2101 piola1(1,2)=stress(6);
2102 piola1(2,1)=stress(6);
2103
2104 switch(fvd.VariableType())
2105 {
2106 case FieldVariableDescriptor::FVT_stress_mises: return piola1.Mises();
2107 case FieldVariableDescriptor::FVT_stress: return fvd.GetComponent(piola1);
2108 case FieldVariableDescriptor::FVT_total_strain: return fvd.GetComponent(strain);
2109 }
2110 }
2111 else // SM-Formulation
2112 {
2113 double xi = p0.X();
2114 Vector3D gamma; GetGamma(xi, gamma, flagD);
2115 Vector3D kappa; GetKappa(xi, kappa, flagD);
2116 Vector3D thicknessstrain; GetThicknessStrain(xi, thicknessstrain, flagD);
2117
2118 //compute interpolated values for gamma, otherwise gamma has oscillatory solution in plot!!!! JG2013-09
2119 //integration rule must be same as in EvalF2SM() !!!!
2120 int order_x_bend = 3; //3
2121 int order_x_thickness_eps = 4; //4
2122 int order_x_axial_eps = 3; //3
2123 int order_x_shear = 3; //3
2124
2125 Vector wT,xT; //integration points and weights
2126 GetIntegrationRule(xT,wT,order_x_axial_eps);
2127
2128 double x1 = GetLx()*0.5*xT(1);
2129 double x2 = GetLx()*0.5*xT(2);
2130 Vector3D gamma_1; GetGamma(x1, gamma_1, flagD);
2131 // Vector3D kappa_1; GetKappa(x1, kappa_1, flagD);
2132 Vector3D gamma_2; GetGamma(x2, gamma_2, flagD);
2133 // Vector3D kappa_2; GetKappa(x2, kappa_2, flagD);
2134
2135 int use_interpolated_values = 1;
2136 //compute interpolated values:
2137 if (use_interpolated_values)
2138 {
2139 gamma = 1./(x2-x1)*((x2-xi)*gamma_1+(xi-x1)*gamma_2); //linear interpolation of all gamma values
2140 }
2141
2142 double beamEA = GetMaterial().BeamEA();
2143 double beamGAky = GetMaterial().BeamGAky();
2144 double beamGAkz = GetMaterial().BeamGAkz();
2145 double beamGAkyz = BeamGAkyz(beamGAky, beamGAkz);
2146 double beamGJkx = GetMaterial().BeamGJkx();
2147 double beamEIy = GetMaterial().BeamEIy();
2148 double beamEIz = GetMaterial().BeamEIz();
2149
2150 Vector3D moment(kappa);
2151 moment.X() *= beamGJkx;
2152 moment.Y() *= beamEIy;
2153 moment.Z() *= beamEIz;
2154
2155 Vector3D beam_force(gamma);
2156 beam_force.X() *= beamEA;
2157 beam_force.Y() *= beamGAky;
2158 beam_force.Z() *= beamGAkz;
2159
2160 /*
2161 //strains and stresses are computed in WRONG WAY JG2013-09 !!!
2162 //the part of bending strain has been forgotten!!!
2163
2164 // for sm-formulation the total strain contains:
2165 // Gamma_1 Gamma_2 Gamma_3
2166 // Gamma_2 E_yy E_yz
2167 // Gamma_3 E_yz E_zz
2168 // ... same for the stress
2169
2170 Matrix3D strains(
2171 gamma.X(), gamma.Y(), gamma.Z(),
2172 gamma.Y(), thicknessstrain.X(), 0.5*thicknessstrain.Z(),
2173 gamma.Z(), 0.5*thicknessstrain.Z(), thicknessstrain.Y()
2174 );
2175
2176 Matrix3D stresses(
2177 beamEA*gamma.X(), beamGAky*gamma.Y(), beamGAkz*gamma.Z(),
2178 beamGAky*gamma.Y(), beamEA*thicknessstrain.X(), beamGAkyz*0.5*thicknessstrain.Z(),
2179 beamGAkz*gamma.Z(), beamGAkyz*0.5*thicknessstrain.Z(), beamEA*thicknessstrain.Y()
2180 );
2181 */
2182
2183 switch(fvd.VariableType())
2184 {
2185 //beam_curvature is defined as 2-component vector! JG2013-09
2186 //case FieldVariableDescriptor::FVT_beam_curvature: return fvd.GetComponent(Vector3D(kappa.X(), kappa.Y(), kappa.Z()));
2187
2188 //strains and stresses are computed in WRONG WAY JG2013-09 !!!
2189 //case FieldVariableDescriptor::FVT_total_strain: return fvd.GetComponent(strains);
2190 //case FieldVariableDescriptor::FVT_stress: return fvd.GetComponent(stresses);
2191
2192 case FieldVariableDescriptor::FVT_beam_torsion: return kappa.X();
2193 case FieldVariableDescriptor::FVT_beam_curvature: return fvd.GetComponent(kappa);
2194
2195 case FieldVariableDescriptor::FVT_beam_moment_torsional: return moment.X();
2196 case FieldVariableDescriptor::FVT_beam_moment_bending: return fvd.GetComponent(moment);
2197
2198 case FieldVariableDescriptor::FVT_beam_axial_extension: return gamma.X();
2199 case FieldVariableDescriptor::FVT_beam_shear: return fvd.GetComponent(gamma);
2200
2201 case FieldVariableDescriptor::FVT_beam_force_axial: return beam_force.X();
2202 case FieldVariableDescriptor::FVT_beam_force_transversal: return fvd.GetComponent(beam_force);
2203 }
2204 }
2205
2206 return FIELD_VARIABLE_NO_VALUE;
2207 }
2208
2209
GetAvailableFieldVariables(TArray<FieldVariableDescriptor> & variables)2210 void ANCFBeamShear3DGeneric::GetAvailableFieldVariables(TArray<FieldVariableDescriptor> & variables)
2211 {
2212 FieldVariableDescriptor::AddTypeIntoArray(variables, FieldVariableDescriptor::FVT_position, FieldVariableDescriptor::FVCI_z);
2213 FieldVariableDescriptor::AddTypeIntoArray(variables, FieldVariableDescriptor::FVT_displacement, FieldVariableDescriptor::FVCI_z);
2214 if (InContinuumMechanicsFormulation())
2215 {
2216 FieldVariableDescriptor::AddTypeIntoArray(variables, FieldVariableDescriptor::FVT_total_strain,
2217 FieldVariableDescriptor::FVCI_z, true);
2218 FieldVariableDescriptor::AddTypeIntoArray(variables, FieldVariableDescriptor::FVT_stress,
2219 FieldVariableDescriptor::FVCI_z, true);
2220 FieldVariableDescriptor::AddTypeIntoArray(variables, FieldVariableDescriptor::FVT_stress_mises);
2221 }
2222 else
2223 {
2224
2225 FieldVariableDescriptor::AddTypeIntoArray(variables, FieldVariableDescriptor::FVT_beam_curvature, FieldVariableDescriptor::FVCI_z);
2226 FieldVariableDescriptor::AddTypeIntoArray(variables, FieldVariableDescriptor::FVT_beam_torsion);
2227 FieldVariableDescriptor::AddTypeIntoArray(variables, FieldVariableDescriptor::FVT_beam_moment_bending, FieldVariableDescriptor::FVCI_z);
2228 FieldVariableDescriptor::AddTypeIntoArray(variables, FieldVariableDescriptor::FVT_beam_moment_torsional);
2229 FieldVariableDescriptor::AddTypeIntoArray(variables, FieldVariableDescriptor::FVT_beam_shear, FieldVariableDescriptor::FVCI_z);
2230 FieldVariableDescriptor::AddTypeIntoArray(variables, FieldVariableDescriptor::FVT_beam_axial_extension);
2231 FieldVariableDescriptor::AddTypeIntoArray(variables, FieldVariableDescriptor::FVT_beam_force_transversal, FieldVariableDescriptor::FVCI_z);
2232 FieldVariableDescriptor::AddTypeIntoArray(variables, FieldVariableDescriptor::FVT_beam_force_axial);
2233
2234 //this is the global beam moment! JG2013-09:
2235 //FieldVariableDescriptor::AddTypeIntoArray(variables, FieldVariableDescriptor::FVT_beam_moment, FieldVariableDescriptor::FVCI_z);
2236 // for sm-formulation the total strain contains:
2237 // Gamma_1 Gamma_2 Gamma_3
2238 // Gamma_2 E_yy E_yz
2239 // Gamma_3 E_yz E_zz
2240 // ... same for the stress
2241
2242 //strains and stresses are computed in WRONG WAY JG2013-09 !!!
2243 //FieldVariableDescriptor::AddTypeIntoArray(variables, FieldVariableDescriptor::FVT_total_strain, FieldVariableDescriptor::FVCI_z, true);
2244 //FieldVariableDescriptor::AddTypeIntoArray(variables, FieldVariableDescriptor::FVT_stress, FieldVariableDescriptor::FVCI_z, true);
2245
2246 }
2247 }
2248
2249 //ACHTUNG: nur kopiert von ANCFBeam3D und ComputeStress-Aufrufe auskommentiert
DrawElement()2250 void ANCFBeamShear3DGeneric::DrawElement()
2251 {
2252 mbs->SetColor(col);
2253
2254 double lx1 = 1; double ly1 = 1/**GetMBS()->GetMagnifyYZ()*/; double lz1 = 1/**GetMBS()->GetMagnifyYZ()*/; //GetMagnifyYZ() leads to wrong stress evaluation
2255 double def_scale = GetMBS()->GetDOption(105); //deformation scaling
2256
2257 int linemode = 1; //0=no lines, 1=outline+color, 2=outline, 3=elementline+color, 4=elementline
2258 if (GetMBS()->GetIOption(110) && !GetMBS()->GetIOption(111))
2259 {
2260 linemode = 2;
2261 }
2262 if (!GetMBS()->GetIOption(110) && GetMBS()->GetIOption(111))
2263 {
2264 linemode = 0;
2265 }
2266 if (!GetMBS()->GetIOption(110) && !GetMBS()->GetIOption(111))
2267 {
2268 linemode = 4;
2269 }
2270
2271 int colormode = 0;
2272 if (GetMBS()->GetActualPostProcessingFieldVariable() != NULL) colormode = 1;
2273 //if (GetMBS()->GetDrawMode() != 0) colormode = 1;//Yuri
2274 else
2275 {
2276 colormode = 0;
2277 }
2278
2279 if (colormode)
2280 {
2281 double modeval = 0;
2282 int xgset = 0;
2283
2284 double tilex = GetMBS()->GetIOption(137);
2285 double tiley = GetMBS()->GetIOption(138);
2286 TArray<Vector3D> points((int)(tilex+1)*(int)(tiley+1));
2287 TArray<double> vals((int)(tilex+1)*(int)(tiley+1));
2288 double v=0;
2289
2290 for (int side=1; side <= 6; side++)
2291 {
2292 points.SetLen(0); vals.SetLen(0);
2293 Vector3D p0, vx, vy;
2294 int tileyn = (int)tiley;
2295 int tilexn = (int)tilex;
2296
2297 if (side == 1)
2298 { //bottom
2299 p0 = Vector3D(-lx1,-ly1,-lz1);
2300 vx = Vector3D(2.*lx1/tilexn,0,0);
2301 vy = Vector3D(0,0,2.*lz1/tileyn);
2302 }
2303 else if (side == 2)
2304 { //top
2305 p0 = Vector3D(-lx1, ly1, lz1);
2306 vx = Vector3D(2.*lx1/tilexn,0,0);
2307 vy = Vector3D(0,0,-2.*lz1/tileyn);
2308 }
2309 else if (side == 3)
2310 { //front
2311 p0 = Vector3D(-lx1,-ly1, lz1);
2312 vx = Vector3D(2.*lx1/tilexn,0,0);
2313 vy = Vector3D(0,2.*ly1/tileyn,0);
2314 }
2315 else if (side == 4)
2316 { //back
2317 p0 = Vector3D(-lx1, ly1,-lz1);
2318 vx = Vector3D(2.*lx1/tilexn,0,0);
2319 vy = Vector3D(0,-2.*ly1/tileyn,0);
2320 }
2321 else if (side == 5)
2322 { //left
2323 p0 = Vector3D(-lx1, ly1,-lz1);
2324 vx = Vector3D(0,-2.*ly1/tilexn,0);
2325 vy = Vector3D(0,0,2.*lz1/tileyn);
2326 }
2327 else if (side == 6)
2328 { //right
2329 p0 = Vector3D( lx1,-ly1,-lz1);
2330 vx = Vector3D(0,2.*ly1/tilexn,0);
2331 vy = Vector3D(0,0,2.*lz1/tileyn);
2332 }
2333
2334 for (double iy = 0; iy <= tileyn+1e-10; iy++)
2335 {
2336 for (double ix = 0; ix <= tilexn+1e-10; ix++)
2337 {
2338 Vector3D ploc = (p0+ix*vx+iy*vy);//-1<=p0<=+1
2339 Vector3D ploc_scaled(ploc(1)*GetLx()*0.5,ploc(2)*GetLy()*0.5,ploc(3)*GetLz()*0.5);//-l/2<=ploc_scaled<=+l/2
2340 Vector3D pg = GetPos3D0D(ploc, def_scale);
2341 //Vector3D p(pg);//p=pg
2342 points.Add(pg);
2343 if (colormode)
2344 v = GetFieldVariableValue(*GetMBS()->GetActualPostProcessingFieldVariable(), ploc_scaled, true);
2345 vals.Add(v);
2346
2347 //if (colormode) ComputeStressD(ploc,type,comp1,comp2,v);//old_KN
2348 //vals.Add(v);
2349 }
2350 }
2351 mbs->DrawColorQuads(points,vals,(int)tilexn+1,(int)tileyn+1,colormode,linemode);
2352 //mbs->UO()<<"points"<<points(1)<<"\n";
2353 }
2354 }
2355 else
2356 {
2357 int drawgrid = GetMBS()->GetIOption(110);
2358 double thickness = GetMBS()->GetDOption(102);
2359 double tiling = GetMBS()->GetIOption(137);
2360
2361 mbs->SetDrawlines(0);
2362 mbs->SetDrawlinesH(0);
2363 Vector3D p1,p2,p3,p4,p5,p6,p7,p8;
2364 for (double i = 0; i < tiling; i++)
2365 {
2366 double l1 = -lx1+2*lx1*i/tiling;
2367 double l2 = -lx1+2*lx1*(i+1)/tiling;
2368 if (i == 0)
2369 {
2370 //p8 = Vector3D(GetPos3D0D(Vector3D(l1,-ly1,-lz1)));
2371 //p7 = Vector3D(GetPos3D0D(Vector3D(l1,-ly1, lz1)));
2372 //p4 = Vector3D(GetPos3D0D(Vector3D(l1, ly1,-lz1)));
2373 //p3 = Vector3D(GetPos3D0D(Vector3D(l1, ly1, lz1)));
2374 p8 = Vector3D(GetPos3D0D(Vector3D(l1,-ly1,-lz1),def_scale));
2375 p7 = Vector3D(GetPos3D0D(Vector3D(l1,-ly1, lz1),def_scale));
2376 p4 = Vector3D(GetPos3D0D(Vector3D(l1, ly1,-lz1),def_scale));
2377 p3 = Vector3D(GetPos3D0D(Vector3D(l1, ly1, lz1),def_scale));
2378 if (drawgrid)
2379 {
2380 mbs->MyDrawLine(p8,p7,thickness);
2381 mbs->MyDrawLine(p7,p3,thickness);
2382 mbs->MyDrawLine(p4,p3,thickness);
2383 mbs->MyDrawLine(p4,p8,thickness);
2384 }
2385 }
2386 else
2387 {
2388 p8 = p6;
2389 p7 = p5;
2390 p4 = p2;
2391 p3 = p1;
2392 }
2393 //p6 = Vector3D(GetPos3D0D(Vector3D(l2,-ly1,-lz1)));
2394 //p5 = Vector3D(GetPos3D0D(Vector3D(l2,-ly1, lz1)));
2395 //p2 = Vector3D(GetPos3D0D(Vector3D(l2, ly1,-lz1)));
2396 //p1 = Vector3D(GetPos3D0D(Vector3D(l2, ly1, lz1)));
2397 p6 = Vector3D(GetPos3D0D(Vector3D(l2,-ly1,-lz1),def_scale));
2398 p5 = Vector3D(GetPos3D0D(Vector3D(l2,-ly1, lz1),def_scale));
2399 p2 = Vector3D(GetPos3D0D(Vector3D(l2, ly1,-lz1),def_scale));
2400 p1 = Vector3D(GetPos3D0D(Vector3D(l2, ly1, lz1),def_scale));
2401 int dout = 0;
2402 if (i == 0 || i == tiling-1) dout = 1;
2403 if (drawgrid) //draw mesh
2404 {
2405 mbs->MyDrawLine(p6,p5,thickness);
2406 mbs->MyDrawLine(p5,p1,thickness);
2407 mbs->MyDrawLine(p2,p1,thickness);
2408 mbs->MyDrawLine(p2,p6,thickness);
2409 mbs->MyDrawLine(p6,p8,thickness);
2410 mbs->MyDrawLine(p5,p7,thickness);
2411 mbs->MyDrawLine(p2,p4,thickness);
2412 mbs->MyDrawLine(p1,p3,thickness);
2413 }
2414 if (GetMBS()->GetIOption(111)) //draw solid
2415 {
2416 mbs->DrawHex(p1,p2,p3,p4,p5,p6,p7,p8,dout);
2417 }
2418 }
2419 mbs->SetDrawlines(0);
2420
2421 //// draw slope vectors
2422 //Vector3D position, slope1, slope2;
2423 //double drawlength = max(ly,lz)*2.;
2424 //double linethickness = max(2.,0.04);
2425 //Vector3D color(0., 0.5, 0.);
2426 //
2427 //position = GetPos3D0D(Vector3D(-1., 0., 0.)); //left node
2428 //slope1(1) = XGD(4) + q0(4);
2429 //slope1(2) = XGD(5) + q0(5);
2430 //slope1(3) = XGD(6) + q0(6);
2431 //slope1.Normalize();
2432 //slope1 *= drawlength;
2433 //mbs->MyDrawLine(position, position+slope1, linethickness, color);
2434 //slope2(1) = XGD(7) + q0(7);
2435 //slope2(2) = XGD(8) + q0(8);
2436 //slope2(3) = XGD(9) + q0(9);
2437 //slope2.Normalize();
2438 //slope2 *= drawlength;
2439 //mbs->MyDrawLine(position, position+slope2, linethickness, color);
2440
2441 //position = GetPos3D0D(Vector3D(1., 0., 0.)); //right node
2442 //slope1(1) = XGD(13) + q0(13);
2443 //slope1(2) = XGD(14) + q0(14);
2444 //slope1(3) = XGD(15) + q0(15);
2445 //slope1.Normalize();
2446 //slope1 *= drawlength;
2447 //mbs->MyDrawLine(position, position+slope1, linethickness, color);
2448 //slope2(1) = XGD(16) + q0(16);
2449 //slope2(2) = XGD(17) + q0(17);
2450 //slope2(3) = XGD(18) + q0(18);
2451 //slope2.Normalize();
2452 //slope2 *= drawlength;
2453 //mbs->MyDrawLine(position, position+slope2, linethickness, color);
2454
2455 //if (nnodes == 3)
2456 //{
2457 // position = GetPos3D0D(Vector3D(0., 0., 0.)); //middle node
2458 // slope1(1) = XGD(22) + q0(22);
2459 // slope1(2) = XGD(23) + q0(23);
2460 // slope1(3) = XGD(24) + q0(24);
2461 // slope1.Normalize();
2462 // slope1 *= drawlength;
2463 // mbs->MyDrawLine(position, position+slope1, linethickness, color);
2464 // slope2(1) = XGD(25) + q0(25);
2465 // slope2(2) = XGD(26) + q0(26);
2466 // slope2(3) = XGD(27) + q0(27);
2467 // slope2.Normalize();
2468 // slope2 *= drawlength;
2469 // mbs->MyDrawLine(position, position+slope2, linethickness, color);
2470 //}
2471 }
2472 };
2473
2474 #pragma endregion
2475
2476 #pragma region ANCFBeamShear3DLinear
2477
ElementDefaultConstructorInitialization()2478 void ANCFBeamShear3DLinear::ElementDefaultConstructorInitialization()
2479 {
2480 // set nodes
2481 nnodes=2;
2482 n1=1; n2=2;
2483
2484 //set reference configuration
2485 q0.SetLen(SOS());
2486 q0.SetAll(0.);
2487 size = Vector3D(1,0.1,0.1);
2488 penalty.SetLen(9);
2489 penalty.SetAll(1);
2490 x_init.SetLen(2*SOS());
2491 x_init.SetAll(0.);
2492
2493 xg.SetLen(SOS());
2494 xg.SetAll(0.);
2495
2496 is_straight_beam_in_reference_configuration = 0;
2497 }
2498
2499 // default set function!!!
SetANCFBeamShear3DLinear(int n1i,int n2i,int materialnumi,const Vector3D & coli)2500 void ANCFBeamShear3DLinear::SetANCFBeamShear3DLinear(int n1i, int n2i, int materialnumi, const Vector3D& coli)
2501 {
2502 ElementDefaultConstructorInitialization();
2503
2504 // set nodes
2505 n1=n1i; n2=n2i;
2506
2507 // set common stuff for ancfbeamshear3d elements
2508 SetANCFBeamShear3DGeneric(materialnumi, coli);
2509
2510 //used for computation speedup: test if beam element is straight in reference configuration
2511 is_straight_beam_in_reference_configuration = 1; // linear (2-noded) element is always straight, only cross section may vary
2512 };
2513
2514 // deprecated set function!!!
SetANCFBeamShear3DLinear(const Vector & xc1,const Vector & xc2,int n1i,int n2i,int materialnumi,const Vector3D & si,const Vector3D & coli)2515 void ANCFBeamShear3DLinear::SetANCFBeamShear3DLinear(const Vector& xc1, const Vector& xc2, int n1i, int n2i, int materialnumi,
2516 const Vector3D& si, const Vector3D& coli)
2517 {
2518 ElementDefaultConstructorInitialization();
2519 // set nodes
2520 //nnodes=2; // is now in ElementDefaultConstructorInitialization
2521 n1=n1i; n2=n2i;
2522
2523 //set reference configuration
2524 //q0.SetLen(SOS()); // is now in ElementDefaultConstructorInitialization
2525 //q0.SetAll(0.); // is now in ElementDefaultConstructorInitialization
2526 int i;
2527 int ndofs_per_node = SOS()/nnodes; //each of the nodes owns ndofs_per_node degrees of freedom
2528 for(int i=1;i<=ndofs_per_node;i++)
2529 {
2530 q0(i) = xc1(i);
2531 q0(i+ndofs_per_node) = xc2(i);
2532 }
2533
2534 // set common stuff for ancfbeamshear3d elements
2535 SetANCFBeamShear3DGeneric(materialnumi, si, coli);
2536
2537 //used for computation speedup: test if beam element is straight in reference configuration
2538 is_straight_beam_in_reference_configuration = 1; // linear (2-noded) element is always straight, only cross section may vary
2539 };
2540
LinkToElements()2541 void ANCFBeamShear3DLinear::LinkToElements()
2542 {
2543 if (SOSowned() == 0)
2544 {
2545 //UO() << "Link nodes to elements in ANCFBeamShear3DLinear\n";
2546 const Node& node1 = GetMBS()->GetNode(n1);
2547 const Node& node2 = GetMBS()->GetNode(n2);
2548 for (int i=1; i <= node1.SOS(); i++)
2549 {
2550 AddLTG(node1.Get(i));
2551 }
2552 for (int i=1; i <= node2.SOS(); i++)
2553 {
2554 AddLTG(node2.Get(i));
2555 }
2556 for (int i=1; i <= node1.SOS(); i++)
2557 {
2558 AddLTG(node1.Get(i+node1.SOS()));
2559 }
2560 for (int i=1; i <= node2.SOS(); i++)
2561 {
2562 AddLTG(node2.Get(i+node2.SOS()));
2563 }
2564 }
2565 }
2566
2567 //-------------------------------------------------------
2568 //Shapefunctions
2569 //
2570 //defined on scaled rectangular element: ploc=(\xi,\eta,\zeta)
2571 //-L/2<\xi<+L/2, -H/2<\eta<+H/2, -B/2<\zeta<+B/2
2572 //-------------------------------------------------------
2573
2574
GetShapes(Vector & sf,const Vector3D & ploc) const2575 void ANCFBeamShear3DLinear::GetShapes(Vector& sf, const Vector3D& ploc) const
2576 //compute all Shapefunctions at ploc=(xi,eta,zeta) in [-L/2,L/2] x [-h/2,h/2] x [-w/2,w/2]
2577 {
2578 sf.SetLen(NS());
2579 double fact = ploc(1)/GetLx();
2580 double sf1 = .5 - fact;
2581 double sf4 = .5 + fact;
2582 sf(1) = sf1;
2583 sf(2) = sf1*ploc(2);
2584 sf(3) = sf1*ploc(3);
2585 sf(4) = sf4;
2586 sf(5) = sf4*ploc(2);
2587 sf(6) = sf4*ploc(3);
2588 }
2589
GetSF(int i,const Vector3D & ploc) const2590 double ANCFBeamShear3DLinear::GetSF(int i, const Vector3D& ploc) const
2591 //computes i-th Shapefunction at ploc=(xi,eta,zeta) in [-L/2,L/2] x [-h/2,h/2] x [-w/2,w/2]
2592 {
2593 switch(i)
2594 {
2595 case 1: return (.5 - ploc(1)/GetLx());
2596 case 2: return (.5 - ploc(1)/GetLx())*ploc(2);
2597 case 3: return (.5 - ploc(1)/GetLx())*ploc(3);
2598 case 4: return (.5 + ploc(1)/GetLx());
2599 case 5: return (.5 + ploc(1)/GetLx())*ploc(2);
2600 case 6: return (.5 + ploc(1)/GetLx())*ploc(3);
2601 default: mbs->UO()<<"only 6 Shapefunctions\n"; return 0.;
2602 }
2603 return 0.;
2604 }
2605
GetShapesx(Vector & sf,const Vector3D & ploc) const2606 void ANCFBeamShear3DLinear::GetShapesx(Vector& sf, const Vector3D& ploc) const
2607 //compute derivative of all Shapefunctions w.r.t. ploc(1)=xi at ploc=(xi,eta,zeta) in [-L/2,L/2] x [-h/2,h/2] x [-w/2,w/2]
2608 {
2609 sf.SetLen(NS());
2610 double fact = 1./GetLx();
2611 sf(1) = -fact;
2612 sf(2) = -fact*ploc(2);
2613 sf(3) = -fact*ploc(3);
2614 sf(4) = fact;
2615 sf(5) = fact*ploc(2);
2616 sf(6) = fact*ploc(3);
2617 }
2618
GetSFx(int i,const Vector3D & ploc) const2619 double ANCFBeamShear3DLinear::GetSFx(int i, const Vector3D& ploc) const
2620 //computes derivative of i-th Shapefunction w.r.t. ploc(1)=xi at ploc=(xi,eta,zeta) in [-L/2,L/2] x [-h/2,h/2] x [-w/2,w/2]
2621 {
2622 switch(i)
2623 {
2624 case 1: return -1./GetLx();
2625 case 2: return -ploc(2)/GetLx();
2626 case 3: return -ploc(3)/GetLx();
2627 case 4: return 1./GetLx();
2628 case 5: return ploc(2)/GetLx();
2629 case 6: return ploc(3)/GetLx();
2630 default: mbs->UO()<<"only 6 Shapefunctions\n"; return 0.;
2631 }
2632 return 0.;
2633 }
2634
GetShapesy(Vector & sf,const Vector3D & ploc) const2635 void ANCFBeamShear3DLinear::GetShapesy(Vector& sf, const Vector3D& ploc) const
2636 //compute derivative of all Shapefunctions w.r.t. ploc(2)=eta at ploc=(xi,eta,zeta) in [-L/2,L/2] x [-h/2,h/2] x [-w/2,w/2]
2637 {
2638 sf.SetLen(NS());
2639 double fact = ploc(1)/GetLx();
2640 sf(1) = 0.;
2641 sf(2) = .5 - fact;
2642 sf(3) = 0.;
2643 sf(4) = 0.;
2644 sf(5) = .5 + fact;
2645 sf(6) = 0.;
2646 }
2647
GetSFy(int i,const Vector3D & ploc) const2648 double ANCFBeamShear3DLinear::GetSFy(int i, const Vector3D& ploc) const
2649 //computes derivative of i-th Shapefunction w.r.t. ploc(2)=eta at ploc=(xi,eta,zeta) in [-L/2,L/2] x [-h/2,h/2] x [-w/2,w/2]
2650 {
2651 switch(i)
2652 {
2653 case 1: return 0.;
2654 case 2: return .5 - ploc(1)/GetLx();
2655 case 3: return 0.;
2656 case 4: return 0.;
2657 case 5: return .5 + ploc(1)/GetLx();
2658 case 6: return 0.;
2659 default: mbs->UO()<<"only 6 Shapefunctions\n"; return 0.;
2660 }
2661 return 0.;
2662 }
2663
GetShapesz(Vector & sf,const Vector3D & ploc) const2664 void ANCFBeamShear3DLinear::GetShapesz(Vector& sf, const Vector3D& ploc) const
2665 //compute derivative of all Shapefunctions w.r.t. ploc(3)=zeta at ploc=(xi,eta,zeta) in [-L/2,L/2] x [-h/2,h/2] x [-w/2,w/2]
2666 {
2667 sf.SetLen(NS());
2668 double fact = ploc(1)/GetLx();
2669 sf(1) = 0.;
2670 sf(2) = 0.;
2671 sf(3) = .5 - fact;
2672 sf(4) = 0.;
2673 sf(5) = 0.;
2674 sf(6) = .5 + fact;
2675 }
2676
GetSFz(int i,const Vector3D & ploc) const2677 double ANCFBeamShear3DLinear::GetSFz(int i, const Vector3D& ploc) const
2678 //computes derivative of i-th Shapefunction w.r.t. ploc(3)=zeta at ploc=(xi,eta,zeta) in [-L/2,L/2] x [-h/2,h/2] x [-w/2,w/2]
2679 {
2680 switch(i)
2681 {
2682 case 1: return 0.;
2683 case 2: return 0.;
2684 case 3: return .5 - ploc(1)/GetLx();
2685 case 4: return 0.;
2686 case 5: return 0.;
2687 case 6: return .5 + ploc(1)/GetLx();
2688 default: mbs->UO()<<"only 6 Shapefunctions\n"; return 0.;
2689 }
2690 return 0.;
2691 }
2692
GetShapesxy(Vector & sf,const Vector3D & ploc) const2693 void ANCFBeamShear3DLinear::GetShapesxy(Vector& sf, const Vector3D& ploc) const
2694 //compute derivative of all Shapefunctions w.r.t. ploc(1)=xi and ploc(2)=eta at ploc=(xi,eta,zeta) in [-L/2,L/2] x [-h/2,h/2] x [-w/2,w/2]
2695 {
2696 sf.SetLen(NS());
2697 double fact = 1./GetLx();
2698 sf(1) = 0.;
2699 sf(2) = -fact;
2700 sf(3) = 0.;
2701 sf(4) = 0.;
2702 sf(5) = fact;
2703 sf(6) = 0.;
2704 }
2705
GetSFxy(int i,const Vector3D & ploc) const2706 double ANCFBeamShear3DLinear::GetSFxy(int i, const Vector3D& ploc) const
2707 //computes derivative of i-th Shapefunction w.r.t. plox(1)=xi and ploc(2)=eta at ploc=(xi,eta,zeta) in [-L/2,L/2] x [-h/2,h/2] x [-w/2,w/2]
2708 {
2709 switch(i)
2710 {
2711 case 1: return 0.;
2712 case 2: return -1./GetLx();
2713 case 3: return 0.;
2714 case 4: return 0.;
2715 case 5: return 1./GetLx();
2716 case 6: return 0.;
2717 default: mbs->UO()<<"only 6 Shapefunctions\n"; return 0.;
2718 }
2719 return 0.;
2720 }
2721
GetShapesxz(Vector & sf,const Vector3D & ploc) const2722 void ANCFBeamShear3DLinear::GetShapesxz(Vector& sf, const Vector3D& ploc) const
2723 //compute derivative of all Shapefunctions w.r.t. ploc(1)=xi and ploc(3)=zeta at ploc=(xi,eta,zeta) in [-L/2,L/2] x [-h/2,h/2] x [-w/2,w/2]
2724 {
2725 sf.SetLen(NS());
2726 double fact = 1./GetLx();
2727 sf(1) = 0.;
2728 sf(2) = 0.;
2729 sf(3) = -fact;
2730 sf(4) = 0.;
2731 sf(5) = 0.;
2732 sf(6) = fact;
2733 }
2734
GetSFxz(int i,const Vector3D & ploc) const2735 double ANCFBeamShear3DLinear::GetSFxz(int i, const Vector3D& ploc) const
2736 //computes derivative of i-th Shapefunction w.r.t. plox(1)=xi and ploc(3)=zeta at ploc=(xi,eta,zeta) in [-L/2,L/2] x [-h/2,h/2] x [-w/2,w/2]
2737 {
2738 switch(i)
2739 {
2740 case 1: return 0.;
2741 case 2: return 0.;
2742 case 3: return -1./GetLx();
2743 case 4: return 0.;
2744 case 5: return 0.;
2745 case 6: return 1./GetLx();
2746 default: mbs->UO()<<"only 6 Shapefunctions\n"; return 0.;
2747 }
2748 return 0.;
2749 }
2750 #pragma endregion
2751
2752
2753
2754 #pragma region ANCFBeamShear3DQuadratic
2755
ElementDefaultConstructorInitialization()2756 void ANCFBeamShear3DQuadratic::ElementDefaultConstructorInitialization()
2757 {
2758 // set nodes
2759 nnodes=3;
2760 n1=1; n2=2; n3 = 3;
2761
2762 //set reference configuration
2763 q0.SetLen(SOS());
2764 q0.SetAll(0.);
2765
2766 //lx = 1;
2767 //ly = 0.1;
2768 //lz = 0.1;
2769 size = Vector3D(1,0.1,0.1);
2770
2771 penalty.SetLen(9);
2772 penalty.SetAll(1);
2773
2774 x_init.SetLen(2*SOS());
2775 x_init.SetAll(0.);
2776
2777 xg.SetLen(SOS());
2778 xg.SetAll(0.);
2779
2780 is_straight_beam_in_reference_configuration = 0;
2781 }
2782
2783 // default set function
SetANCFBeamShear3DQuadratic(int n1i,int n2i,int n3i,int materialnumi,const Vector3D & coli)2784 void ANCFBeamShear3DQuadratic::SetANCFBeamShear3DQuadratic(int n1i, int n2i, int n3i, int materialnumi, const Vector3D& coli)
2785 {
2786 ElementDefaultConstructorInitialization();
2787 // set nodes
2788 //nnodes=3; // is now in ElementDefaultConstructorInitialization()
2789 n1=n1i; n2=n2i, n3=n3i;
2790
2791 // set common stuff for ancfbeamshear3d elements
2792 SetANCFBeamShear3DGeneric(materialnumi, coli);
2793
2794 //used for computation speedup: test if beam element is straight in reference configuration
2795 Vector3D pos_left(q0(1), q0(2), q0(3));
2796 Vector3D pos_right(q0(10), q0(11), q0(12));
2797 Vector3D pos_mid(q0(19), q0(20), q0(21));
2798 Vector3D v1 = pos_right - pos_mid;
2799 Vector3D v2 = pos_mid - pos_left;
2800 v1.Normalize();
2801 v2.Normalize();
2802 is_straight_beam_in_reference_configuration = (v1*v2 > 1-1e-14); // if both sections of the beam are colinear, then this scalar product gives numerically 1
2803 };
2804
2805 // deprecated set function!!!
SetANCFBeamShear3DQuadratic(const Vector & xc1,const Vector & xc2,const Vector & xc3,int n1i,int n2i,int n3i,int materialnumi,const Vector3D & si,const Vector3D & coli)2806 void ANCFBeamShear3DQuadratic::SetANCFBeamShear3DQuadratic(const Vector& xc1, const Vector& xc2, const Vector& xc3, int n1i, int n2i, int n3i, int materialnumi,
2807 const Vector3D& si, const Vector3D& coli)
2808 {
2809 ElementDefaultConstructorInitialization();
2810 // set nodes
2811 //nnodes=3; // is now in ElementDefaultConstructorInitialization()
2812 n1=n1i; n2=n2i, n3=n3i;
2813
2814 //set reference configuration
2815 //q0.SetLen(SOS()); // is now in ElementDefaultConstructorInitialization()
2816 //q0.SetAll(0.); // is now in ElementDefaultConstructorInitialization()
2817 int i;
2818 int ndofs_per_node = SOS()/nnodes; //each of the nodes owns ndofs_per_node degrees of freedom
2819 for(int i=1;i<=ndofs_per_node;i++)
2820 {
2821 q0(i) = xc1(i);
2822 q0(i+ndofs_per_node) = xc2(i);
2823 q0(i+2*ndofs_per_node) = xc3(i);
2824 }
2825
2826 //used for computation speedup: test if beam element is straight in reference configuration
2827 Vector3D pos_left(q0(1), q0(2), q0(3));
2828 Vector3D pos_right(q0(10), q0(11), q0(12));
2829 Vector3D pos_mid(q0(10), q0(11), q0(12));
2830 Vector3D v1 = pos_right - pos_mid;
2831 Vector3D v2 = pos_mid - pos_left;
2832 v1.Normalize();
2833 v2.Normalize();
2834 is_straight_beam_in_reference_configuration = (v1*v2 > 1-1e-14); // if both sections of the beam are colinear, then this scalar product gives numerically 1
2835
2836 // set common stuff for ancfbeamshear3d elements
2837 SetANCFBeamShear3DGeneric(materialnumi, si, coli);
2838 };
2839
LinkToElements()2840 void ANCFBeamShear3DQuadratic::LinkToElements()
2841 {
2842 if (SOSowned() == 0)
2843 {
2844 //UO() << "Link nodes to elements in ANCFBeamShear3DQuadratic\n";
2845 const Node& node1 = GetMBS()->GetNode(n1);
2846 const Node& node2 = GetMBS()->GetNode(n2);
2847 const Node& node3 = GetMBS()->GetNode(n3);
2848 for (int i=1; i <= node1.SOS(); i++)
2849 {
2850 AddLTG(node1.Get(i));
2851 }
2852 for (int i=1; i <= node2.SOS(); i++)
2853 {
2854 AddLTG(node2.Get(i));
2855 }
2856 for (int i=1; i <= node3.SOS(); i++)
2857 {
2858 AddLTG(node3.Get(i));
2859 }
2860 for (int i=1; i <= node1.SOS(); i++)
2861 {
2862 AddLTG(node1.Get(i+node1.SOS()));
2863 }
2864 for (int i=1; i <= node2.SOS(); i++)
2865 {
2866 AddLTG(node2.Get(i+node2.SOS()));
2867 }
2868 for (int i=1; i <= node3.SOS(); i++)
2869 {
2870 AddLTG(node3.Get(i+node3.SOS()));
2871 }
2872 }
2873 }
2874
2875 //-------------------------------------------------------
2876 //Shapefunctions
2877 //
2878 //defined on scaled rectangular element: ploc=(\xi,\eta,\zeta)
2879 //-L/2<\xi<+L/2, -H/2<\eta<+H/2, -B/2<\zeta<+B/2
2880 //-------------------------------------------------------
2881
2882
GetShapes(Vector & sf,const Vector3D & ploc) const2883 void ANCFBeamShear3DQuadratic::GetShapes(Vector& sf, const Vector3D& ploc) const
2884 //compute all Shapefunctions at ploc=(xi,eta,zeta) in [-L/2,L/2] x [-h/2,h/2] x [-w/2,w/2]
2885 {
2886 sf.SetLen(NS()); //NS=3*nnodes, vector sf has 9 entries
2887 double fact=2./(GetLx()*GetLx());
2888 sf(1)=-fact*(GetLx()/2.-ploc(1))*ploc(1);
2889 sf(2)=sf(1)*ploc(2);
2890 sf(3)=sf(1)*ploc(3);
2891 sf(4)=fact*(GetLx()/2.+ploc(1))*ploc(1);
2892 sf(5)=sf(4)*ploc(2);
2893 sf(6)=sf(4)*ploc(3);
2894 sf(7)=-2.*fact*(GetLx()/2.+ploc(1))*(ploc(1)-GetLx()/2.);
2895 sf(8)=sf(7)*ploc(2);
2896 sf(9)=sf(7)*ploc(3);
2897 }
2898
GetSF(int i,const Vector3D & ploc) const2899 double ANCFBeamShear3DQuadratic::GetSF(int i, const Vector3D& ploc) const
2900 //computes i-th Shapefunction at ploc=(xi,eta,zeta) in [-L/2,L/2] x [-h/2,h/2] x [-w/2,w/2]
2901 {
2902 switch(i)
2903 {
2904 case 1: return -2./(GetLx()*GetLx())*ploc(1)*(GetLx()/2.-ploc(1));
2905 case 2: return -2./(GetLx()*GetLx())*ploc(1)*ploc(2)*(GetLx()/2.-ploc(1));
2906 case 3: return -2./(GetLx()*GetLx())*ploc(1)*ploc(3)*(GetLx()/2.-ploc(1));
2907 case 4: return 2./(GetLx()*GetLx())*ploc(1)*(GetLx()/2.+ploc(1));
2908 case 5: return 2./(GetLx()*GetLx())*ploc(1)*ploc(2)*(GetLx()/2.+ploc(1));
2909 case 6: return 2./(GetLx()*GetLx())*ploc(1)*ploc(3)*(GetLx()/2.+ploc(1));
2910 case 7: return -4./(GetLx()*GetLx())*(ploc(1)-GetLx()/2.)*(ploc(1)+GetLx()/2.);
2911 case 8: return -4./(GetLx()*GetLx())*ploc(2)*(ploc(1)-GetLx()/2.)*(ploc(1)+GetLx()/2.);
2912 case 9: return -4./(GetLx()*GetLx())*ploc(3)*(ploc(1)-GetLx()/2.)*(ploc(1)+GetLx()/2.);
2913 default: mbs->UO()<<"only 9 Shapefunctions\n"; return 0.;
2914 }
2915 return 0.;
2916 }
2917
GetShapesx(Vector & sf,const Vector3D & ploc) const2918 void ANCFBeamShear3DQuadratic::GetShapesx(Vector& sf, const Vector3D& ploc) const
2919 //compute derivative of all Shapefunctions w.r.t. ploc(1)=xi at ploc=(xi,eta,zeta) in [-L/2,L/2] x [-h/2,h/2] x [-w/2,w/2]
2920 {
2921 sf.SetLen(NS()); //NS=3*nnodes, vector sf has 9 entries
2922 double fact = 4.*ploc(1)/(GetLx()*GetLx());
2923 double sf1 = -1./GetLx()+fact;
2924 double sf4 = 1./GetLx()+fact;
2925 double sf7 = -2.*fact;
2926 sf(1) = sf1;
2927 sf(2) = sf1*ploc(2);
2928 sf(3) = sf1*ploc(3);
2929 sf(4) = sf4;
2930 sf(5) = sf4*ploc(2);
2931 sf(6) = sf4*ploc(3);
2932 sf(7) = sf7;
2933 sf(8) = sf7*ploc(2);
2934 sf(9) = sf7*ploc(3);
2935 }
2936
GetSFx(int i,const Vector3D & ploc) const2937 double ANCFBeamShear3DQuadratic::GetSFx(int i, const Vector3D& ploc) const
2938 //computes derivative of i-th Shapefunction w.r.t. ploc(1)=xi at ploc=(xi,eta,zeta) in [-L/2,L/2] x [-h/2,h/2] x [-w/2,w/2]
2939 {
2940 switch(i)
2941 {
2942 case 1: return -1./GetLx()+4.*ploc(1)/(GetLx()*GetLx());
2943 case 2: return -ploc(2)/GetLx()+4.*ploc(1)*ploc(2)/(GetLx()*GetLx());
2944 case 3: return -ploc(3)/GetLx()+4.*ploc(1)*ploc(3)/(GetLx()*GetLx());
2945 case 4: return 1./GetLx()+4.*ploc(1)/(GetLx()*GetLx());
2946 case 5: return ploc(2)/GetLx()+4.*ploc(1)*ploc(2)/(GetLx()*GetLx());
2947 case 6: return ploc(3)/GetLx()+4.*ploc(1)*ploc(3)/(GetLx()*GetLx());
2948 case 7: return -8.*ploc(1)/(GetLx()*GetLx());
2949 case 8: return -8.*ploc(1)*ploc(2)/(GetLx()*GetLx());
2950 case 9: return -8.*ploc(1)*ploc(3)/(GetLx()*GetLx());
2951 default: mbs->UO()<<"only 9 Shapefunctions\n"; return 0.;
2952 }
2953 return 0.;
2954 }
2955
GetShapesy(Vector & sf,const Vector3D & ploc) const2956 void ANCFBeamShear3DQuadratic::GetShapesy(Vector& sf, const Vector3D& ploc) const
2957 //compute derivative of all Shapefunctions w.r.t. ploc(2)=eta at ploc=(xi,eta,zeta) in [-L/2,L/2] x [-h/2,h/2] x [-w/2,w/2]
2958 {
2959 sf.SetLen(NS()); //NS=3*nnodes, vector sf has 9 entries
2960 double fact1 = ploc(1)/GetLx();
2961 double fact2 = 2.*fact1*fact1;
2962 sf(1) = 0;
2963 sf(2) = -fact1 + fact2;
2964 sf(3) = 0;
2965 sf(4) = 0;
2966 sf(5) = fact1 + fact2;
2967 sf(6) = 0;
2968 sf(7) = 0;
2969 sf(8) = 1. - 2.*fact2;
2970 sf(9) = 0;
2971 }
2972
GetSFy(int i,const Vector3D & ploc) const2973 double ANCFBeamShear3DQuadratic::GetSFy(int i, const Vector3D& ploc) const
2974 //computes derivative of i-th Shapefunction w.r.t. ploc(2)=eta at ploc=(xi,eta,zeta) in [-L/2,L/2] x [-h/2,h/2] x [-w/2,w/2]
2975 {
2976 switch(i)
2977 {
2978 case 1: return 0;
2979 case 2: return -ploc(1)/GetLx()+2.*ploc(1)*ploc(1)/(GetLx()*GetLx());
2980 case 3: return 0;
2981 case 4: return 0;
2982 case 5: return ploc(1)/GetLx()+2.*ploc(1)*ploc(1)/(GetLx()*GetLx());
2983 case 6: return 0;
2984 case 7: return 0;
2985 case 8: return 1.-4.*ploc(1)*ploc(1)/(GetLx()*GetLx());
2986 case 9: return 0;
2987 default: mbs->UO()<<"only 9 Shapefunctions\n"; return 0.;
2988 }
2989 return 0.;
2990 }
2991
GetShapesxy(Vector & sf,const Vector3D & ploc) const2992 void ANCFBeamShear3DQuadratic::GetShapesxy(Vector& sf, const Vector3D& ploc) const
2993 //compute derivative of all Shapefunctions w.r.t. ploc(1)=xi and ploc(2)=eta at ploc=(xi,eta,zeta) in [-L/2,L/2] x [-h/2,h/2] x [-w/2,w/2]
2994 {
2995 sf.SetLen(NS()); //NS=3*nnodes, vector sf has 9 entries
2996 double fact1 = 1./GetLx();
2997 double fact2 = 4.*ploc(1)*fact1*fact1;
2998 sf(1) = 0;
2999 sf(2) = -fact1 + fact2;
3000 sf(3) = 0;
3001 sf(4) = 0;
3002 sf(5) = fact1 + fact2;
3003 sf(6) = 0;
3004 sf(7) = 0;
3005 sf(8) = -2.*fact2;
3006 sf(9) = 0;
3007 }
3008
GetSFxy(int i,const Vector3D & ploc) const3009 double ANCFBeamShear3DQuadratic::GetSFxy(int i, const Vector3D& ploc) const
3010 //computes derivative of i-th Shapefunction w.r.t. ploc(1)=xi and ploc(2)=eta at ploc=(xi,eta,zeta) in [-L/2,L/2] x [-h/2,h/2] x [-w/2,w/2]
3011 {
3012 switch(i)
3013 {
3014 case 1: return 0;
3015 case 2: return -1./GetLx()+4.*ploc(1)/(GetLx()*GetLx());
3016 case 3: return 0;
3017 case 4: return 0;
3018 case 5: return 1./GetLx()+4.*ploc(1)/(GetLx()*GetLx());
3019 case 6: return 0;
3020 case 7: return 0;
3021 case 8: return -8.*ploc(1)/(GetLx()*GetLx());
3022 case 9: return 0;
3023 default: mbs->UO()<<"only 9 Shapefunctions\n"; return 0.;
3024 }
3025 return 0.;
3026 }
3027
GetShapesz(Vector & sf,const Vector3D & ploc) const3028 void ANCFBeamShear3DQuadratic::GetShapesz(Vector& sf, const Vector3D& ploc) const
3029 //compute derivative of all Shapefunctions w.r.t. ploc(3)=zeta at ploc=(xi,eta,zeta) in [-L/2,L/2] x [-h/2,h/2] x [-w/2,w/2]
3030 {
3031 sf.SetLen(NS()); //NS=3*nnodes, vector sf has 9 entries
3032 double fact1 = ploc(1)/GetLx();
3033 double fact2 = 2.*fact1*fact1;
3034 sf(1) = 0;
3035 sf(2) = 0;
3036 sf(3) = -fact1 + fact2;
3037 sf(4) = 0;
3038 sf(5) = 0;
3039 sf(6) = fact1 + fact2;
3040 sf(7) = 0;
3041 sf(8) = 0;
3042 sf(9) = 1. - 2.*fact2;
3043 }
3044
GetSFz(int i,const Vector3D & ploc) const3045 double ANCFBeamShear3DQuadratic::GetSFz(int i, const Vector3D& ploc) const
3046 //computes derivative of i-th Shapefunction w.r.t. ploc(3)=zeta at ploc=(xi,eta,zeta) in [-L/2,L/2] x [-h/2,h/2] x [-w/2,w/2]
3047 {
3048 switch(i)
3049 {
3050 case 1: return 0;
3051 case 2: return 0;
3052 case 3: return -ploc(1)/GetLx()+2.*ploc(1)*ploc(1)/(GetLx()*GetLx());
3053 case 4: return 0;
3054 case 5: return 0;
3055 case 6: return ploc(1)/GetLx()+2.*ploc(1)*ploc(1)/(GetLx()*GetLx());
3056 case 7: return 0;
3057 case 8: return 0;
3058 case 9: return 1.-4.*ploc(1)*ploc(1)/(GetLx()*GetLx());
3059 default: mbs->UO()<<"only 9 Shapefunctions\n"; return 0.;
3060 }
3061 return 0.;
3062 }
3063
GetShapesxz(Vector & sf,const Vector3D & ploc) const3064 void ANCFBeamShear3DQuadratic::GetShapesxz(Vector& sf, const Vector3D& ploc) const
3065 //compute derivative of all Shapefunctions w.r.t. ploc(1)=xi and ploc(3)=zeta at ploc=(xi,eta,zeta) in [-L/2,L/2] x [-h/2,h/2] x [-w/2,w/2]
3066 {
3067 sf.SetLen(NS()); //NS=3*nnodes, vector sf has 9 entries
3068 double fact1 = 1./GetLx();
3069 double fact2 = 4.*ploc(1)*fact1*fact1;
3070 sf(1) = 0;
3071 sf(2) = 0;
3072 sf(3) = -fact1 + fact2;
3073 sf(4) = 0;
3074 sf(5) = 0;
3075 sf(6) = fact1 + fact2;
3076 sf(7) = 0;
3077 sf(8) = 0;
3078 sf(9) = -2.*fact2;
3079 }
3080
GetSFxz(int i,const Vector3D & ploc) const3081 double ANCFBeamShear3DQuadratic::GetSFxz(int i, const Vector3D& ploc) const
3082 //computes derivative of i-th Shapefunction w.r.t. ploc(1)=xi and ploc(2)=zeta at ploc=(xi,eta,zeta) in [-L/2,L/2] x [-h/2,h/2] x [-w/2,w/2]
3083 {
3084 switch(i)
3085 {
3086 case 1: return 0;
3087 case 2: return 0;
3088 case 3: return -1./GetLx()+4.*ploc(1)/(GetLx()*GetLx());
3089 case 4: return 0;
3090 case 5: return 0;
3091 case 6: return 1./GetLx()+4.*ploc(1)/(GetLx()*GetLx());
3092 case 7: return 0;
3093 case 8: return 0;
3094 case 9: return -8.*ploc(1)/(GetLx()*GetLx());
3095 default: mbs->UO()<<"only 9 Shapefunctions\n"; return 0.;
3096 }
3097 return 0.;
3098 }
3099
3100 #pragma endregion
3101