beam_elements.cc
Go to the documentation of this file.
1// LIC// ====================================================================
2// LIC// This file forms part of oomph-lib, the object-oriented,
3// LIC// multi-physics finite-element library, available
4// LIC// at http://www.oomph-lib.org.
5// LIC//
6// LIC// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
7// LIC//
8// LIC// This library is free software; you can redistribute it and/or
9// LIC// modify it under the terms of the GNU Lesser General Public
10// LIC// License as published by the Free Software Foundation; either
11// LIC// version 2.1 of the License, or (at your option) any later version.
12// LIC//
13// LIC// This library is distributed in the hope that it will be useful,
14// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16// LIC// Lesser General Public License for more details.
17// LIC//
18// LIC// You should have received a copy of the GNU Lesser General Public
19// LIC// License along with this library; if not, write to the Free Software
20// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21// LIC// 02110-1301 USA.
22// LIC//
23// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24// LIC//
25// LIC//====================================================================
26// Non-inline functions for Kirchhoff Love beam elements
27
28#include "beam_elements.h"
29
30namespace oomph
31{
32 //================================================================
33 /// Static default value for 2nd Piola Kirchhoff prestress (zero)
34 //================================================================
36
37 //=====================================================================
38 /// Static default value for timescale ratio (1.0 for natural scaling)
39 //=====================================================================
41
42 //=========================================================
43 /// Static default value for non-dim wall thickness (1/20)
44 //=========================================================
45 // i.e. the reference thickness 'h_0'
47
48 //=======================================================================
49 /// Default load function (zero traction)
50 //=======================================================================
52 const Vector<double>& x,
53 const Vector<double>& N,
54 Vector<double>& load)
55 {
56 unsigned n_dim = load.size();
57 for (unsigned i = 0; i < n_dim; i++)
58 {
59 load[i] = 0.0;
60 }
61 }
62
63 //=======================================================================
64 /// Default wall profile function (constant thickness h_0)
65 //=======================================================================
67 const Vector<double>& x,
68 double& h_ratio)
69 {
70 h_ratio = 1.0;
71 }
72
73
74 //=======================================================================
75 /// Get position vector to and normal vector on wall
76 //=======================================================================
80 {
81#ifdef PARANOID
82 if (N.size() != 2)
83 {
84 std::ostringstream error_message;
85 error_message << "Normal vector should have dimension 2, not" << N.size()
86 << std::endl;
87
88 throw OomphLibError(
89 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
90 }
91 if (r.size() != 2)
92 {
93 std::ostringstream error_message;
94 error_message << "Position vector should have dimension 2, not"
95 << r.size() << std::endl;
96
97 throw OomphLibError(
98 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
99 }
100
101 if (s.size() != 1)
102 {
103 std::ostringstream error_message;
104 error_message << "Local coordinate should have dimension 1, not"
105 << s.size() << std::endl;
106
107 throw OomphLibError(
108 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
109 }
110#endif
111
112 // Set the dimension of the global coordinates
113 unsigned n_dim = Undeformed_beam_pt->ndim();
114
115 // Set the number of lagrangian coordinates
116 unsigned n_lagrangian = Undeformed_beam_pt->nlagrangian();
117
118 // Find out how many nodes there are
119 unsigned n_node = nnode();
120
121 // Find out how many positional dofs there are
122 unsigned n_position_type = nnodal_position_type();
123
124 // Set up memory for the shape functions:
125
126 // # of nodes, # of positional dofs
127 Shape psi(n_node, n_position_type);
128
129 // # of nodes, # of positional dofs, # of lagrangian coords (for deriv)
130 DShape dpsidxi(n_node, n_position_type, n_lagrangian);
131
132 // Call the derivatives of the shape functions w.r.t. Lagrangian coords
133 dshape_lagrangian(s, psi, dpsidxi);
134
135 // Base Vector
136 DenseMatrix<double> interpolated_A(n_lagrangian, n_dim);
137
138 // Initialise to zero
139
140 // Loop over coordinate directions/components of Vector
141 for (unsigned i = 0; i < n_dim; i++)
142 {
143 r[i] = 0.0;
144 // Loop over derivatives/base Vectors (just one here)
145 for (unsigned j = 0; j < n_lagrangian; j++)
146 {
147 interpolated_A(j, i) = 0.0;
148 }
149 }
150
151 // Loop over directions
152 for (unsigned i = 0; i < n_dim; i++)
153 {
154 // Loop over nodes
155 for (unsigned l = 0; l < n_node; l++)
156 {
157 // Loop over positional dofs
158 for (unsigned k = 0; k < n_position_type; k++)
159 {
160 r[i] += raw_nodal_position_gen(l, k, i) * psi(l, k);
161
162 // Loop over derivative directions (just one here)
163 for (unsigned j = 0; j < n_lagrangian; j++)
164 {
165 interpolated_A(j, i) +=
166 raw_nodal_position_gen(l, k, i) * dpsidxi(l, k, j);
167 }
168 }
169 }
170 }
171
172 // Calculate the length of the normal vector
173 double length = pow(interpolated_A(0, 0), 2) + pow(interpolated_A(0, 1), 2);
174
175 // Calculate the normal
176 N[0] = -interpolated_A(0, 1) / sqrt(length);
177 N[1] = interpolated_A(0, 0) / sqrt(length);
178 }
179
180
181 //=======================================================================
182 /// Get position vector to and non-unit tangent vector on wall: dr/ds
183 //=======================================================================
186 Vector<double>& drds)
187 {
188#ifdef PARANOID
189 if (drds.size() != 2)
190 {
191 std::ostringstream error_message;
192 error_message << "Tangent vector should have dimension 2, not"
193 << drds.size() << std::endl;
194
195 throw OomphLibError(
196 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
197 }
198 if (r.size() != 2)
199 {
200 std::ostringstream error_message;
201 error_message << "Position vector should have dimension 2, not"
202 << r.size() << std::endl;
203
204 throw OomphLibError(
205 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
206 }
207
208 if (s.size() != 1)
209 {
210 std::ostringstream error_message;
211 error_message << "Local coordinate should have dimension 1, not"
212 << s.size() << std::endl;
213
214 throw OomphLibError(
215 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
216 }
217#endif
218
219 // Set the dimension of the global coordinates
220 unsigned n_dim = Undeformed_beam_pt->ndim();
221
222 // Set the number of lagrangian coordinates
223 unsigned n_lagrangian = Undeformed_beam_pt->nlagrangian();
224
225 // Find out how many nodes there are
226 unsigned n_node = nnode();
227
228 // Find out how many positional dofs there are
229 unsigned n_position_type = nnodal_position_type();
230
231 // Set up memory for the shape functions:
232
233 // # of nodes, # of positional dofs
234 Shape psi(n_node, n_position_type);
235
236 // # of nodes, # of positional dofs, # of lagrangian coords (for deriv)
237 DShape dpsids(n_node, n_position_type, n_lagrangian);
238
239 // Call the derivatives of the shape functions w.r.t. local coords
240 dshape_local(s, psi, dpsids);
241
242 // Initialise to zero
243
244 // Loop over coordinate directions/components of Vector
245 for (unsigned i = 0; i < n_dim; i++)
246 {
247 r[i] = 0.0;
248 drds[i] = 0.0;
249 }
250
251 // Loop over directions
252 for (unsigned i = 0; i < n_dim; i++)
253 {
254 // Loop over nodes
255 for (unsigned l = 0; l < n_node; l++)
256 {
257 // Loop over positional dofs
258 for (unsigned k = 0; k < n_position_type; k++)
259 {
260 r[i] += raw_nodal_position_gen(l, k, i) * psi(l, k);
261 // deriv w.r.t. to zero-th (and only) local coordiate
262 drds[i] += raw_nodal_position_gen(l, k, i) * dpsids(l, k, 0);
263 }
264 }
265 }
266 }
267
268
269 //=======================================================================
270 /// Return the residuals for the equations of Kirchhoff-Love beam
271 /// theory with linear constitutive equations; if Solid_ic_pt!=0, we
272 /// assign residuals which force the assignement of an initial shape/
273 /// veloc/accel to the dofs.
274 //=======================================================================
276 Vector<double>& residuals)
277 {
278 // Set up the initial conditions, if an IC pointer has been set
279 if (Solid_ic_pt != 0)
280 {
282 return;
283 }
284
285 // Set the dimension of the global coordinates
286 const unsigned n_dim = Undeformed_beam_pt->ndim();
287
288 // Set the number of lagrangian coordinates
289 const unsigned n_lagrangian = Undeformed_beam_pt->nlagrangian();
290
291 // Find out how many nodes there are
292 const unsigned n_node = nnode();
293
294 // Find out how many positional dofs there are
295 const unsigned n_position_type = nnodal_position_type();
296
297 // Integer to store the local equation number
298 int local_eqn = 0;
299
300 // Setup memory for accelerations
301 Vector<double> accel(2);
302
303 // Set up memory for the shape functions:
304
305 // # of nodes, # of positional dofs
306 Shape psi(n_node, n_position_type);
307
308 // # of nodes, # of positional dofs, # of lagrangian coords (for deriv)
309 DShape dpsidxi(n_node, n_position_type, n_lagrangian);
310
311 // # of nodes, # of positional dofs, # of derivs)
312 DShape d2psidxi(n_node, n_position_type, n_lagrangian);
313
314 // Set # of integration points
315 const unsigned n_intpt = integral_pt()->nweight();
316
317 // Get Physical Variables from Element
318
319 // Thickness h_0/R, axial prestress, timescale ratio (density)
320 const double HoR_0 = h(); // i.e. refers to reference thickness 'h_0'
321 const double sigma_0 = sigma0();
322 const double Lambda_sq = lambda_sq();
323
324 // Loop over the integration points
325 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
326 {
327 // Get the integral weight
328 double w = integral_pt()->weight(ipt);
329
330 // Call the derivatives of the shape functions w.r.t. Lagrangian coords
331 double J = d2shape_lagrangian_at_knot(ipt, psi, dpsidxi, d2psidxi);
332
333 // Premultiply the weights and the Jacobian
334 double W = w * J;
335
336 // Calculate local values of lagrangian position and
337 // the derivative of global position wrt lagrangian coordinates
338 Vector<double> interpolated_xi(n_lagrangian, 0.0),
339 interpolated_x(n_dim, 0.0);
340 DenseMatrix<double> interpolated_A(n_lagrangian, n_dim);
341 DenseMatrix<double> interpolated_dAdxi(n_lagrangian, n_dim);
342
343 // Initialise to zero
344 // Loop over coordinate directions/components of Vector
345 for (unsigned i = 0; i < n_dim; i++)
346 {
347 // Initialise acclerations
348 accel[i] = 0.0;
349 // Loop over derivatives/base Vectors (just one here)
350 for (unsigned j = 0; j < n_lagrangian; j++)
351 {
352 interpolated_A(j, i) = 0.0;
353 }
354 // Loop over derivatives of base Vector (just one here)
355 for (unsigned j = 0; j < n_lagrangian; j++)
356 {
357 interpolated_dAdxi(j, i) = 0.0;
358 }
359 }
360
361 // Calculate displacements, accelerations and spatial derivatives
362 for (unsigned l = 0; l < n_node; l++)
363 {
364 // Loop over positional dofs
365 for (unsigned k = 0; k < n_position_type; k++)
366 {
367 // Loop over Lagrangian coordinate directions [xi_gen[] are the
368 // the *gen*eralised Lagrangian coordinates: node, type, direction]
369 for (unsigned i = 0; i < n_lagrangian; i++)
370 {
372 raw_lagrangian_position_gen(l, k, i) * psi(l, k);
373 }
374
375 // Loop over components of the deformed position Vector
376 for (unsigned i = 0; i < n_dim; i++)
377 {
378 interpolated_x[i] += raw_nodal_position_gen(l, k, i) * psi(l, k);
379
380 accel[i] += raw_dnodal_position_gen_dt(2, l, k, i) * psi(l, k);
381
382 // Loop over derivative directions (just one here)
383 for (unsigned j = 0; j < n_lagrangian; j++)
384 {
385 interpolated_A(j, i) +=
386 raw_nodal_position_gen(l, k, i) * dpsidxi(l, k, j);
387 }
388
389 // Loop over the second derivative directions (just one here)
390 for (unsigned j = 0; j < n_lagrangian; j++)
391 {
392 interpolated_dAdxi(j, i) +=
393 raw_nodal_position_gen(l, k, i) * d2psidxi(l, k, j);
394 }
395 }
396 }
397 }
398
399 // Setup position Vector and derivatives of undeformed config
400 Vector<double> R(n_dim);
401 DenseMatrix<double> a(n_lagrangian, n_dim);
402 RankThreeTensor<double> dadxi(n_lagrangian, n_lagrangian, n_dim);
403
404 // Get the undeformed geometry
406
407 // Declare and calculate the undeformed and deformed metric tensor
408 // and the strain tensor (these are just scalars)
409 double amet = 0.0, Amet = 0.0;
410
411 // Now calculate the dot product
412 for (unsigned k = 0; k < n_dim; k++)
413 {
414 amet += a(0, k) * a(0, k);
415 Amet += interpolated_A(0, k) * interpolated_A(0, k);
416 }
417
418 double gamma = 0.5 * (Amet - amet);
419
420 // Calculate the contravariant metric tensors
421 double adet = amet; // double aup = 1.0/amet;
422 double Adet = Amet; // double Aup = 1.0/Amet;
423
424 // Calculate the unit normal Vectors
425 Vector<double> n(2);
426 Vector<double> N(2);
427 n[0] = -a(0, 1) / sqrt(adet);
428 n[1] = a(0, 0) / sqrt(adet);
429
430 N[0] = -interpolated_A(0, 1) / sqrt(Adet);
431 N[1] = interpolated_A(0, 0) / sqrt(Adet);
432
433 // Square root of deformed determinant
434 double sqrt_Adet = sqrt(Adet);
435
436 // Calculate the curvature tensors
437 double b = n[0] * dadxi(0, 0, 0) + n[1] * dadxi(0, 0, 1);
438 double B =
439 N[0] * interpolated_dAdxi(0, 0) + N[1] * interpolated_dAdxi(0, 1);
440
441 // Set up the change of curvature tensor
442 double kappa = b - B;
443
444 // Define the load Vector
445 Vector<double> f(n_dim);
446
447 // Get load Vector
449
450 // Define the wall thickness ratio profile
451 double h_ratio = 0;
452
453 // Get wall thickness ratio profile
455
456 // Thickness h/R
457 double HoR = HoR_0 * h_ratio;
458
459 // Allocate storage for variations of normal Vector
460 double normal_var[n_dim][n_dim];
461
462 // Geometrically nonlinear beam equations from
463 //--------------------------------------------
464 // principle of virtual displacements
465 //-----------------------------------
466
467 // Loop over the number of nodes
468 for (unsigned n = 0; n < n_node; n++)
469 {
470 // Loop over the type of degree of freedom
471 for (unsigned k = 0; k < n_position_type; k++)
472 {
473 // Setup components of variation of normal Vector:
474 // ...(variation w.r.t. direction, component)
475
476 normal_var[0][0] = interpolated_A(0, 0) * interpolated_A(0, 1) *
477 dpsidxi(n, k, 0) /
478 (sqrt_Adet * sqrt_Adet * sqrt_Adet);
479
480 normal_var[1][0] =
481 (interpolated_A(0, 1) * interpolated_A(0, 1) / Adet - 1.0) *
482 dpsidxi(n, k, 0) / (sqrt_Adet);
483
484 normal_var[0][1] =
485 (1.0 - interpolated_A(0, 0) * interpolated_A(0, 0) / Adet) *
486 dpsidxi(n, k, 0) / (sqrt_Adet);
487
488 normal_var[1][1] = -interpolated_A(0, 0) * interpolated_A(0, 1) *
489 dpsidxi(n, k, 0) /
490 (sqrt_Adet * sqrt_Adet * sqrt_Adet);
491
492 // Loop over the coordinate directions
493 for (unsigned i = 0; i < n_dim; i++)
494 {
495 // Find the equation number
496 local_eqn = position_local_eqn(n, k, i);
497
498 // If it's not a boundary condition
499 if (local_eqn >= 0)
500 {
501 // Premultiply by thickness profile
502
503 // External forcing
504 residuals[local_eqn] -=
505 h_ratio * (1.0 / HoR) * f[i] * psi(n, k) * W * sqrt(Adet);
506
507 residuals[local_eqn] +=
508 h_ratio * Lambda_sq * accel[i] * psi(n, k) * W * sqrt(adet);
509
510 // Membrane term with axial prestress
511 residuals[local_eqn] += h_ratio * (sigma_0 + gamma) *
512 interpolated_A(0, i) * dpsidxi(n, k, 0) *
513 W * sqrt(adet);
514
515 // Bending term: Minus sign because \delta \kappa = - \delta B
516 residuals[local_eqn] -=
517 h_ratio * (1.0 / 12.0) * HoR * HoR * kappa *
518 (N[i] * d2psidxi(n, k, 0) +
519 normal_var[i][0] * interpolated_dAdxi(0, 0) +
520 normal_var[i][1] * interpolated_dAdxi(0, 1)) *
521 W * sqrt(adet);
522 }
523 }
524 }
525 }
526
527 } // End of loop over the integration points
528 }
529
530
531 //=======================================================================
532 /// Get FE jacobian and residuals (Jacobian done by finite differences)
533 //=======================================================================
535 Vector<double>& residuals, DenseMatrix<double>& jacobian)
536 {
537 // Call the element's residuals vector
539
540 // Solve for the consistent acceleration in Newmark scheme?
542 {
544 return;
545 }
546
547 // Allocate storage for the full residuals of the element
548 unsigned n_dof = ndof();
549 Vector<double> full_residuals(n_dof);
550
551 // Call the full residuals
552 get_residuals(full_residuals);
553
554 // Get the solid entries in the jacobian using finite differences
556 full_residuals, jacobian);
557
558 // Get the entries from the external data, usually load terms
560 jacobian);
561 }
562
563
564 //=======================================================================
565 /// Compute the potential (strain) and kinetic energy of the
566 /// element (wrapper function).
567 //=======================================================================
568 void KirchhoffLoveBeamEquations::get_energy(double& pot_en, double& kin_en)
569 {
570 // Initialise
571 double stretch = 0.0;
572 double bend = 0.0;
573 kin_en = 0.0;
574
575 // Compute energy
576 get_energy(stretch, bend, kin_en);
577
578 // Sum components of potential energy
579 pot_en = stretch + bend;
580 }
581
582
583 //=======================================================================
584 /// Compute the potential (strain) and kinetic energy of the
585 /// element, breaking down the potential energy into stretching and
586 /// bending components.
587 //=======================================================================
589 double& bend,
590 double& kin_en)
591 {
592 // Initialise
593 stretch = 0.0;
594 bend = 0.0;
595 kin_en = 0.0;
596
597 // Set the dimension of the global coordinates
598 unsigned n_dim = Undeformed_beam_pt->ndim();
599
600 // Set the number of lagrangian coordinates
601 unsigned n_lagrangian = Undeformed_beam_pt->nlagrangian();
602
603 // Find out how many nodes there are
604 unsigned n_node = nnode();
605
606 // Find out how many positional dofs there are
607 unsigned n_position_dofs = nnodal_position_type();
608
609 // Setup memory for veloc
610 Vector<double> veloc(2);
611
612 // Set up memory for the shape functions:
613
614 // # of nodes, # of positional dofs
615 Shape psi(n_node, n_position_dofs);
616
617 // # of nodes, # of positional dofs, # of lagrangian coords (for deriv)
618 DShape dpsidxi(n_node, n_position_dofs, n_lagrangian);
619
620 // # of nodes, # of positional dofs, # of derivs)
621 DShape d2psidxi(n_node, n_position_dofs, n_lagrangian);
622
623 // Set # of integration points
624 unsigned n_intpt = integral_pt()->nweight();
625
626 // Get Physical Variables from Element
627
628 // Thickness h_0/R, axial prestress, timescale ratio (density)
629 double HoR_0 = h(); // i.e. refers to reference thickness 'h_0'
630 double sigma_0 = sigma0();
631
632 double Lambda_sq = lambda_sq();
633
634 // Loop over the integration points
635 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
636 {
637 // Get the integral weight
638 double w = integral_pt()->weight(ipt);
639
640 // Call the derivatives of the shape functions w.r.t. Lagrangian coords
641 double J = d2shape_lagrangian_at_knot(ipt, psi, dpsidxi, d2psidxi);
642
643 // Premultiply the weights and the Jacobian
644 double W = w * J;
645
646 // Calculate local values of lagrangian position and
647 // the derivative of global position wrt lagrangian coordinates
648 Vector<double> interpolated_xi(n_lagrangian, 0.0),
649 interpolated_x(n_dim, 0.0);
650
651 DenseMatrix<double> interpolated_A(n_lagrangian, n_dim);
652 DenseMatrix<double> interpolated_dAdxi(n_lagrangian, n_dim);
653
654 // Initialise to zero
655
656 // Loop over coordinate directions/components of Vector
657 for (unsigned i = 0; i < n_dim; i++)
658 {
659 // Initialise veloc
660 veloc[i] = 0.0;
661
662 // Loop over derivatives/base Vectors (just one here)
663 for (unsigned j = 0; j < n_lagrangian; j++)
664 {
665 interpolated_A(j, i) = 0.0;
666 }
667 // Loop over derivatives of base Vector (just one here)
668 for (unsigned j = 0; j < n_lagrangian; j++)
669 {
670 interpolated_dAdxi(j, i) = 0.0;
671 }
672 }
673
674 // Calculate displacements, accelerations and spatial derivatives
675 for (unsigned l = 0; l < n_node; l++)
676 {
677 // Loop over positional dofs
678 for (unsigned k = 0; k < n_position_dofs; k++)
679 {
680 // Loop over Lagrangian coordinate directions [xi_gen[] are the
681 // the *gen*eralised Lagrangian coordinates: node, type, direction]
682 for (unsigned i = 0; i < n_lagrangian; i++)
683 {
685 raw_lagrangian_position_gen(l, k, i) * psi(l, k);
686 }
687
688 // Loop over components of the deformed position Vector
689 for (unsigned i = 0; i < n_dim; i++)
690 {
691 // Need this for wall thickness ratio profile
692 interpolated_x[i] += raw_nodal_position_gen(l, k, i) * psi(l, k);
693
694 veloc[i] += raw_dnodal_position_gen_dt(1, l, k, i) * psi(l, k);
695
696 // Loop over derivative directions (just one here)
697 for (unsigned j = 0; j < n_lagrangian; j++)
698 {
699 interpolated_A(j, i) +=
700 raw_nodal_position_gen(l, k, i) * dpsidxi(l, k, j);
701 }
702
703 // Loop over the second derivative directions (just one here)
704 for (unsigned j = 0; j < n_lagrangian; j++)
705 {
706 interpolated_dAdxi(j, i) +=
707 raw_nodal_position_gen(l, k, i) * d2psidxi(l, k, j);
708 }
709 }
710 }
711 }
712
713 // Get square of veloc
714 double veloc_sq = 0;
715 for (unsigned i = 0; i < n_dim; i++)
716 {
717 veloc_sq += veloc[i] * veloc[i];
718 }
719
720 // Setup position Vector and derivatives of undeformed config
721 Vector<double> R(n_dim);
722 DenseMatrix<double> a(n_lagrangian, n_dim);
723 RankThreeTensor<double> dadxi(n_lagrangian, n_lagrangian, n_dim);
724
725 // Get the undeformed geometry
727
728 // Declare and calculate the undeformed and deformed metric tensor
729 // and the strain tensor (these are 1d tensors, i.e. scalars)
730 double amet = 0.0, Amet = 0.0;
731
732 // Work out metric and strain tensors
733 // Now calculate the dot product
734 for (unsigned k = 0; k < n_dim; k++)
735 {
736 amet += a(0, k) * a(0, k);
737 Amet += interpolated_A(0, k) * interpolated_A(0, k);
738 }
739
740 // Calculate strain tensor
741 double gamma = 0.5 * (Amet - amet);
742
743 // Calculate the contravariant metric tensors
744 double adet = amet; // aup = 1.0/amet;
745 double Adet = Amet; // Aup = 1.0/Amet;
746
747 // Calculate the unit normal Vectors
748 Vector<double> n(2);
749 Vector<double> N(2);
750 n[0] = -a(0, 1) / sqrt(adet);
751 n[1] = a(0, 0) / sqrt(adet);
752
753 N[0] = -interpolated_A(0, 1) / sqrt(Adet);
754 N[1] = interpolated_A(0, 0) / sqrt(Adet);
755
756 // Calculate the curvature tensors
757 double b = n[0] * dadxi(0, 0, 0) + n[1] * dadxi(0, 0, 1);
758 double B =
759 N[0] * interpolated_dAdxi(0, 0) + N[1] * interpolated_dAdxi(0, 1);
760
761 // Set up the change of curvature tensor
762 double kappa = b - B;
763
764 // Define the wall thickness ratio profile
765 double h_ratio = 0;
766
767 // Get wall thickness ratio profile
769
770 // Thickness h/R
771 double HoR = HoR_0 * h_ratio;
772
773 // Add contributions
774 stretch +=
775 h_ratio * 0.5 * (gamma + sigma_0) * (gamma + sigma_0) * W * sqrt(adet);
776 bend += h_ratio * 0.5 * (1.0 / 12.0) * HoR * HoR * kappa * kappa * W *
777 sqrt(adet);
778 kin_en += h_ratio * 0.5 * Lambda_sq * veloc_sq * W * sqrt(adet);
779 } // End of loop over the integration points
780 }
781
782 //=======================================================================
783 /// Output position at previous time (t=0: present; t>0: previous)
784 /// with specified number of plot points.
785 //=======================================================================
786 void HermiteBeamElement::output(const unsigned& t,
787 std::ostream& outfile,
788 const unsigned& n_plot) const
789 {
790#ifdef WARN_ABOUT_SUBTLY_CHANGED_OOMPH_INTERFACES
791 // Warn about time argument being moved to the front
793 "Order of function arguments has changed between versions 0.8 and 0.85",
794 "HermiteBeamElement::output(...)",
795 OOMPH_EXCEPTION_LOCATION);
796#endif
797
798 // Local variables
799 Vector<double> s(1);
800
801 // Tecplot header info
802 outfile << "ZONE I=" << n_plot << std::endl;
803
804 // Find the dimension of the first node
805 unsigned n_dim = this->node_pt(0)->ndim();
806
807 // Loop over plot points
808 for (unsigned l = 0; l < n_plot; l++)
809 {
810 s[0] = -1.0 + l * 2.0 / (n_plot - 1);
811
812 // Output the Eulerian coordinates
813 for (unsigned i = 0; i < n_dim; i++)
814 {
815 outfile << interpolated_x(t, s, i) << " ";
816 }
817 outfile << std::endl;
818 }
819 }
820
821
822 //=======================================================================
823 /// Output function
824 //=======================================================================
825 void HermiteBeamElement::output(std::ostream& outfile, const unsigned& n_plot)
826 {
827 // Local variables
828 Vector<double> s(1);
829
830 // Tecplot header info
831 outfile << "ZONE I=" << n_plot << std::endl;
832
833 // Set the number of lagrangian coordinates
834 unsigned n_lagrangian = Undeformed_beam_pt->nlagrangian();
835
836 // Set the dimension of the global coordinates
837 unsigned n_dim = Undeformed_beam_pt->ndim();
838
839 // Find out how many nodes there are
840 unsigned n_node = nnode();
841
842 // Find out how many positional dofs there are
843 unsigned n_position_dofs = nnodal_position_type();
844
845 Vector<double> veloc(n_dim);
846 Vector<double> accel(n_dim);
847 Vector<double> posn(n_dim);
848
849 // # of nodes, # of positional dofs
850 Shape psi(n_node, n_position_dofs);
851
852 // Loop over element plot points
853 for (unsigned l1 = 0; l1 < n_plot; l1++)
854 {
855 s[0] = -1.0 + l1 * 2.0 / (n_plot - 1);
856
857 // Get shape functions
858 shape(s, psi);
859
860 Vector<double> interpolated_xi(n_lagrangian);
861 interpolated_xi[0] = 0.0;
862
863 // Loop over coordinate directions/components of Vector
864 for (unsigned i = 0; i < n_dim; i++)
865 {
866 // Initialise acclerations and veloc
867 accel[i] = 0.0;
868 veloc[i] = 0.0;
869 posn[i] = 0.0;
870 }
871
872
873 // Calculate displacements, accelerations and spatial derivatives
874 for (unsigned l = 0; l < n_node; l++)
875 {
876 // Loop over positional dofs
877 for (unsigned k = 0; k < n_position_dofs; k++)
878 {
879 // Loop over Lagrangian coordinate directions [xi_gen[] are the
880 // the *gen*eralised Lagrangian coordinates: node, type, direction]
881 for (unsigned i = 0; i < n_lagrangian; i++)
882 {
884 raw_lagrangian_position_gen(l, k, i) * psi(l, k);
885 }
886
887 // Loop over components of the deformed position Vector
888 for (unsigned i = 0; i < n_dim; i++)
889 {
890 accel[i] += raw_dnodal_position_gen_dt(2, l, k, i) * psi(l, k);
891 veloc[i] += raw_dnodal_position_gen_dt(1, l, k, i) * psi(l, k);
892 posn[i] += raw_dnodal_position_gen_dt(0, l, k, i) * psi(l, k);
893 }
894 }
895 }
896
897 double scalar_accel = 0.0;
898 double scalar_veloc = 0.0;
899 double scalar_posn = 0.0;
900
901 // Output position etc.
902 for (unsigned i = 0; i < n_dim; i++)
903 {
904 outfile << posn[i] << " ";
905 scalar_posn += pow(posn[i], 2);
906 }
907 outfile << interpolated_xi[0] << " ";
908 for (unsigned i = 0; i < n_dim; i++)
909 {
910 outfile << veloc[i] << " ";
911 scalar_veloc += pow(veloc[i], 2);
912 }
913 for (unsigned i = 0; i < n_dim; i++)
914 {
915 outfile << accel[i] << " ";
916 scalar_accel += pow(accel[i], 2);
917 }
918 outfile << sqrt(scalar_posn) << " ";
919 outfile << sqrt(scalar_veloc) << " ";
920 outfile << sqrt(scalar_accel) << " ";
921 outfile << std::endl;
922 }
923 }
924
925
926 //=======================================================================
927 /// Output function
928 //=======================================================================
929 void HermiteBeamElement::output(std::ostream& outfile)
930 {
931 unsigned n_plot = 5;
932 output(outfile, n_plot);
933 }
934
935
936 //=======================================================================
937 /// C-style output position at previous time (t=0: present; t>0: previous)
938 /// with specified number of plot points.
939 //=======================================================================
940 void HermiteBeamElement::output(const unsigned& t,
941 FILE* file_pt,
942 const unsigned& n_plot) const
943 {
944#ifdef WARN_ABOUT_SUBTLY_CHANGED_OOMPH_INTERFACES
945 // Warn about time argument being moved to the front
947 "Order of function arguments has changed between versions 0.8 and 0.85",
948 "HermiteBeamElement::output(...)",
949 OOMPH_EXCEPTION_LOCATION);
950#endif
951
952 // Local variables
953 Vector<double> s(1);
954
955 // Tecplot header info
956 // outfile << "ZONE I=" << n_plot << std::endl;
957 fprintf(file_pt, "ZONE I=%i \n", n_plot);
958
959 // Find the dimension of the first node
960 unsigned n_dim = this->node_pt(0)->ndim();
961
962 // Loop over plot points
963 for (unsigned l = 0; l < n_plot; l++)
964 {
965 s[0] = -1.0 + l * 2.0 / (n_plot - 1);
966
967 // Output the Eulerian coordinates
968 for (unsigned i = 0; i < n_dim; i++)
969 {
970 // outfile << interpolated_x(s,t,i) << " " ;
971 fprintf(file_pt, "%g ", interpolated_x(t, s, i));
972 }
973 // outfile << std::endl;
974 fprintf(file_pt, "\n");
975 }
976 }
977
978
979 //=======================================================================
980 /// Output function
981 //=======================================================================
982 void HermiteBeamElement::output(FILE* file_pt, const unsigned& n_plot)
983 {
984 // Local variables
985 Vector<double> s(1);
986
987 // Tecplot header info
988 // outfile << "ZONE I=" << n_plot << std::endl;
989 fprintf(file_pt, "ZONE I=%i \n", n_plot);
990
991 // Set the number of lagrangian coordinates
992 unsigned n_lagrangian = Undeformed_beam_pt->nlagrangian();
993
994 // Set the dimension of the global coordinates
995 unsigned n_dim = Undeformed_beam_pt->ndim();
996
997 // Find out how many nodes there are
998 unsigned n_node = nnode();
999
1000 // Find out how many positional dofs there are
1001 unsigned n_position_dofs = nnodal_position_type();
1002
1003 Vector<double> veloc(n_dim);
1004 Vector<double> accel(n_dim);
1005 Vector<double> posn(n_dim);
1006
1007 // # of nodes, # of positional dofs
1008 Shape psi(n_node, n_position_dofs);
1009
1010 // Loop over element plot points
1011 for (unsigned l1 = 0; l1 < n_plot; l1++)
1012 {
1013 s[0] = -1.0 + l1 * 2.0 / (n_plot - 1);
1014
1015 // Get shape functions
1016 shape(s, psi);
1017
1018 Vector<double> interpolated_xi(n_lagrangian);
1019 interpolated_xi[0] = 0.0;
1020
1021 // Loop over coordinate directions/components of Vector
1022 for (unsigned i = 0; i < n_dim; i++)
1023 {
1024 // Initialise acclerations and veloc
1025 accel[i] = 0.0;
1026 veloc[i] = 0.0;
1027 posn[i] = 0.0;
1028 }
1029
1030
1031 // Calculate displacements, accelerations and spatial derivatives
1032 for (unsigned l = 0; l < n_node; l++)
1033 {
1034 // Loop over positional dofs
1035 for (unsigned k = 0; k < n_position_dofs; k++)
1036 {
1037 // Loop over Lagrangian coordinate directions [xi_gen[] are the
1038 // the *gen*eralised Lagrangian coordinates: node, type, direction]
1039 for (unsigned i = 0; i < n_lagrangian; i++)
1040 {
1041 interpolated_xi[i] +=
1042 raw_lagrangian_position_gen(l, k, i) * psi(l, k);
1043 }
1044
1045 // Loop over components of the deformed position Vector
1046 for (unsigned i = 0; i < n_dim; i++)
1047 {
1048 accel[i] += raw_dnodal_position_gen_dt(2, l, k, i) * psi(l, k);
1049 veloc[i] += raw_dnodal_position_gen_dt(1, l, k, i) * psi(l, k);
1050 posn[i] += raw_dnodal_position_gen_dt(0, l, k, i) * psi(l, k);
1051 }
1052 }
1053 }
1054
1055 double scalar_accel = 0.0;
1056 double scalar_veloc = 0.0;
1057 double scalar_posn = 0.0;
1058
1059 // Output position etc.
1060 for (unsigned i = 0; i < n_dim; i++)
1061 {
1062 // outfile << posn[i] << " " ;
1063 fprintf(file_pt, "%g ", posn[i]);
1064 scalar_posn += pow(posn[i], 2);
1065 }
1066 // outfile << interpolated_xi[0] << " ";
1067 fprintf(file_pt, "%g ", interpolated_xi[0]);
1068 for (unsigned i = 0; i < n_dim; i++)
1069 {
1070 // outfile << veloc[i] << " " ;
1071 fprintf(file_pt, "%g ", veloc[i]);
1072 scalar_veloc += pow(veloc[i], 2);
1073 }
1074 for (unsigned i = 0; i < n_dim; i++)
1075 {
1076 // outfile << accel[i] << " " ;
1077 fprintf(file_pt, "%g ", accel[i]);
1078 scalar_accel += pow(accel[i], 2);
1079 }
1080 // outfile << sqrt(scalar_posn) << " ";
1081 fprintf(file_pt, "%g ", sqrt(scalar_posn));
1082 // outfile << sqrt(scalar_veloc) << " ";
1083 fprintf(file_pt, "%g ", sqrt(scalar_veloc));
1084 // outfile << sqrt(scalar_accel) << " ";
1085 fprintf(file_pt, "%g ", sqrt(scalar_accel));
1086 // outfile << std::endl;
1087 fprintf(file_pt, "\n");
1088 }
1089 }
1090
1091
1092 //=======================================================================
1093 /// C-style output function
1094 //=======================================================================
1095 void HermiteBeamElement::output(FILE* file_pt)
1096 {
1097 unsigned n_plot = 5;
1098 output(file_pt, n_plot);
1099 }
1100
1101 /// ////////////////////////////////////////////////////////////////////////
1102 /// ////////////////////////////////////////////////////////////////////////
1103 /// ////////////////////////////////////////////////////////////////////////
1104
1105
1106 //=======================================================================
1107 /// Function used to find the local coordinate s that corresponds to the
1108 /// "global" intrinsic coordinate zeta (in the element's incarnation
1109 /// as a GeomObject). For this element, zeta is equal to the Lagrangian
1110 /// coordinate xi. If the required zeta is located within this
1111 /// element, geom_object_pt points to "this" element. If zeta
1112 /// is not located within the element it is set to NULL.
1113 //======================================================================
1115 const Vector<double>& zeta,
1116 GeomObject*& geom_object_pt,
1118 const bool& use_coordinate_as_initial_guess)
1119 {
1120 // Assumed that the first node has a lower xi coordinate than the second
1121 unsigned lo = 0, hi = 1;
1122
1123 // If the first node has a higher xi then swap
1125 {
1126 lo = 1;
1127 hi = 0;
1128 }
1129
1130 // Tolerance for finding zeta
1131 double epsilon = 1.0e-13;
1132
1133 // If zeta is not in the element, then return a null pointer
1134 // Correct for rounding errors here
1135 if ((zeta[0] - raw_lagrangian_position(lo, 0) < -epsilon) ||
1136 (zeta[0] - raw_lagrangian_position(hi, 0) > epsilon))
1137 {
1138 geom_object_pt = 0;
1139 return;
1140 }
1141
1142 // Otherwise, zeta is located in this element. For now assume
1143 // that the relationship between zeta and s is linear as it will
1144 // be for uniform node spacing... We'll trap this further down.
1145 // In general we'll need a Newton method here and we'll implement
1146 // this as soon as we have an example with non-uniformly spaced
1147 // FSIHermiteBeamElements...
1148
1149 // The GeomObject that contains the zeta coordinate is "this":
1150 geom_object_pt = this;
1151
1152
1153 // Find the fraction along the element
1154 double zeta_fraction =
1155 (zeta[0] - raw_lagrangian_position(lo, 0)) /
1157
1158 s[0] = -1.0 + zeta_fraction * 2.0;
1159
1160#ifdef PARANOID
1161 // Check if we've really ended where we wanted to be...
1162 Vector<double> zeta_test(1);
1163 interpolated_zeta(s, zeta_test);
1164 if (std::fabs(zeta[0] - zeta_test[0]) > epsilon)
1165 {
1166 std::ostringstream error_stream;
1167 error_stream
1168 << "The zeta coordinate " << zeta_test[0] << " \n"
1169 << "computed by interpolated_zeta() for s[0]=" << s[0] << " \n"
1170 << "differs by more than the tolerance (" << epsilon << ") from \n "
1171 << "the required value " << zeta_test[0] << " \n\n"
1172 << "You're probably using a mesh with non-uniformly \n "
1173 << "spaced FSIHermiteBeamElements. For such cases the root finding"
1174 << "in FSIHermiteBeamElement::locate_zeta() must be replaced "
1175 << "by a proper Newton method or some such thing...\n";
1176 OomphLibError(error_stream.str(),
1177 "FSIHermiteBeamElement::locate_zeta()",
1178 OOMPH_EXCEPTION_LOCATION);
1179 }
1180#endif
1181
1182
1183 // Check for rounding error
1184 if (s[0] > 1.0)
1185 {
1186 s[0] = 1.0;
1187 }
1188 if (s[0] < -1.0)
1189 {
1190 s[0] = -1.0;
1191 }
1192 }
1193
1194
1195 //=========================================================================
1196 /// Define the dposition function. This is used to set no-slip boundary
1197 /// conditions in FSI problems.
1198 //=========================================================================
1200 const Vector<double>& s, DenseMatrix<double>& drdxi) const
1201 {
1202#ifdef PARANOID
1203 if (Undeformed_beam_pt == 0)
1204 {
1205 throw OomphLibError(
1206 "Undeformed_beam_pt has not been set",
1207 "FSIHermiteBeamElement::dposition_dlagrangian_at_local_coordinate()",
1208 OOMPH_EXCEPTION_LOCATION);
1209 }
1210#endif
1211
1212 // Find the dimension of the global coordinate
1213 unsigned n_dim = Undeformed_beam_pt->ndim();
1214 // Find the number of lagrangian coordinates
1215 unsigned n_lagrangian = Undeformed_beam_pt->nlagrangian();
1216 // Find out how many nodes there are
1217 unsigned n_node = nnode();
1218 // Find out how many positional dofs there are
1219 unsigned n_position_type = this->nnodal_position_type();
1220
1221 // Set up memory for the shape functions
1222 Shape psi(n_node, n_position_type);
1223 DShape dpsidxi(n_node, n_position_type, n_lagrangian);
1224
1225 // Get the derivatives of the shape functions w.r.t local coordinates
1226 // at this point
1227 dshape_lagrangian(s, psi, dpsidxi);
1228
1229 // Initialise the derivatives to zero
1230 drdxi.initialise(0.0);
1231
1232 // Loop over the nodes
1233 for (unsigned l = 0; l < n_node; l++)
1234 {
1235 // Loop over the positional types
1236 for (unsigned k = 0; k < n_position_type; k++)
1237 {
1238 // Loop over the dimensions
1239 for (unsigned i = 0; i < n_dim; i++)
1240 {
1241 // Loop over the lagrangian coordinates
1242 for (unsigned j = 0; j < n_lagrangian; j++)
1243 {
1244 // Add the contribution to the derivative
1245 drdxi(j, i) += nodal_position_gen(l, k, i) * dpsidxi(l, k, j);
1246 }
1247 }
1248 }
1249 }
1250 }
1251
1252 //=============================================================================
1253 /// Create a list of pairs for all unknowns in this element,
1254 /// so that the first entry in each pair contains the global equation
1255 /// number of the unknown, while the second one contains the number
1256 /// of the "DOF type" that this unknown is associated with.
1257 /// (Function can obviously only be called if the equation numbering
1258 /// scheme has been set up.)
1259 /// This element is only in charge of the solid dofs.
1260 //=============================================================================
1262 std::list<std::pair<unsigned long, unsigned>>& dof_lookup_list) const
1263 {
1264 // temporary pair (used to store dof lookup prior to being added to list)
1265 std::pair<unsigned, unsigned> dof_lookup;
1266
1267 // number of nodes
1268 const unsigned n_node = this->nnode();
1269
1270 // Get the number of position dofs and dimensions at the node
1271 const unsigned n_position_type = nnodal_position_type();
1272 const unsigned nodal_dim = nodal_dimension();
1273
1274 // Integer storage for local unknown
1275 int local_unknown = 0;
1276
1277 // Loop over the nodes
1278 for (unsigned n = 0; n < n_node; n++)
1279 {
1280 // Loop over position dofs
1281 for (unsigned k = 0; k < n_position_type; k++)
1282 {
1283 // Loop over dimension
1284 for (unsigned i = 0; i < nodal_dim; i++)
1285 {
1286 // If the variable is free
1287 local_unknown = position_local_eqn(n, k, i);
1288
1289 // ignore pinned values
1290 if (local_unknown >= 0)
1291 {
1292 // store dof lookup in temporary pair: First entry in pair
1293 // is global equation number; second entry is dof type
1294 dof_lookup.first = this->eqn_number(local_unknown);
1295 dof_lookup.second = 0;
1296
1297 // add to list
1298 dof_lookup_list.push_front(dof_lookup);
1299 }
1300 }
1301 }
1302 }
1303 }
1304
1305
1306 /// ///////////////////////////////////////////////////////////////////
1307 /// ///////////////////////////////////////////////////////////////////
1308 /// ///////////////////////////////////////////////////////////////////
1309
1310
1311 //===========================================================================
1312 /// Constructor, takes the pointer to the "bulk" element, the
1313 /// index of the fixed local coordinate and its value represented
1314 /// by an integer (+/- 1), indicating that the face is located
1315 /// at the max. or min. value of the "fixed" local coordinate
1316 /// in the bulk element.
1317 //===========================================================================
1320 FiniteElement* const& bulk_el_pt, const int& face_index)
1322 {
1323 // Number of nodal position types: 2
1325
1326 // Let the bulk element build the FaceElement, i.e. setup the pointers
1327 // to its nodes (by referring to the appropriate nodes in the bulk
1328 // element), etc.
1329 bulk_el_pt->build_face_element(face_index, this);
1330
1331 // Resize vector to some point on the symmetry line along which the
1332 // end of the beam is sliding. Initialise for symmetry line = y-axis
1333 Vector_to_symmetry_line.resize(2);
1334 Vector_to_symmetry_line[0] = 0.0;
1335 Vector_to_symmetry_line[1] = 0.0;
1336
1337 // Resize normal vector to the symmetry line along which the
1338 // end of the beam is sliding. Initialise for symmetry line = y-axis
1339 Normal_to_symmetry_line.resize(2);
1340 Normal_to_symmetry_line[0] = -1.0;
1341 Normal_to_symmetry_line[1] = 0.0;
1342
1343 // Resize number of dofs at the element's one and only node by two
1344 // to accomodate the Lagrange multipliers
1345
1346 // We need two additional values at the one-and-only node in this element
1347 Vector<unsigned> nadditional_data_values(1);
1348 nadditional_data_values[0] = 2;
1349 resize_nodes(nadditional_data_values);
1350 }
1351
1352
1353 //===========================================================================
1354 /// Add the element's contribution to its residual vector
1355 //===========================================================================
1358 {
1359 // Get position vector to and normal & tangent vectors on wall (from
1360 // bulk element)
1361 HermiteBeamElement* bulk_el_pt =
1362 dynamic_cast<HermiteBeamElement*>(bulk_element_pt());
1363
1364 // Get the value of the constant local coordinate in the bulk element
1365 // Dummy local coordinate of size zero
1366 Vector<double> s(0);
1367 Vector<double> s_bulk(1);
1368 this->get_local_coordinate_in_bulk(s, s_bulk);
1369
1370 // Get position vector to and normal on wall
1371 Vector<double> r(2);
1372 Vector<double> n(2);
1373 bulk_el_pt->get_normal(s_bulk, r, n);
1374
1375 // Get (non-unit) tangent vector
1376 Vector<double> drds(2);
1377 bulk_el_pt->get_non_unit_tangent(s_bulk, r, drds);
1378
1379 // Residual corresponding to: Point must be on symmetry line
1380 double res0 =
1383
1384 // Work out tangent to symmetry line
1385 Vector<double> tangent_to_symmetry_line(2);
1386 tangent_to_symmetry_line[0] = -Normal_to_symmetry_line[1];
1387 tangent_to_symmetry_line[1] = Normal_to_symmetry_line[0];
1388
1389 // Residual corresponding to: Beam must meet symmetry line at right angle
1390 double res1 = drds[0] * tangent_to_symmetry_line[0] +
1391 drds[1] * tangent_to_symmetry_line[1];
1392
1393 // The first Lagrange multiplier enforces the along-the-beam displacement:
1394 int j_local = nodal_local_eqn(0, Nbulk_value[0]);
1395 residuals[j_local] += res0;
1396
1397 // Second Lagrange multiplier enforces the derivative of the
1398 // transverse displacement:
1399 j_local = nodal_local_eqn(0, Nbulk_value[0] + 1);
1400 residuals[j_local] += res1;
1401
1402 // Add Lagrange multiplier contribution to the bulk equations
1403
1404 // Lagrange multipliers
1405 double lambda0 = node_pt(0)->value(Nbulk_value[0]);
1406 double lambda1 = node_pt(0)->value(Nbulk_value[0] + 1);
1407
1408 // Loop over the solid dofs
1409
1410 // How many positional values are there?
1411 unsigned n_dim = nodal_dimension();
1412 unsigned n_type = nnodal_position_type();
1413
1414 /// Loop over directions
1415 for (unsigned i = 0; i < n_dim; i++)
1416 {
1417 // Loop over types
1418 for (unsigned k = 0; k < n_type; k++)
1419 {
1420 // Get local equation number
1421 int j_local = position_local_eqn(0, k, i);
1422
1423 // Real dof?
1424 if (j_local >= 0)
1425 {
1426 // Derivative of first residual w.r.t. to discrete displacements:
1427 // Only the type-zero dofs make a contribution:
1428 double dres0dxik = 0.0;
1429 if (k == 0)
1430 {
1431 dres0dxik = Normal_to_symmetry_line[i];
1432 }
1433
1434 // Derivative of second residual w.r.t. to discrete displacements:
1435 // Only the type-one dofs make a contribution:
1436 double dres1dxik = 0.0;
1437 if (k == 1)
1438 {
1439 dres1dxik = tangent_to_symmetry_line[i];
1440 }
1441
1442 // Add to the residual
1443 residuals[j_local] += lambda0 * dres0dxik + lambda1 * dres1dxik;
1444 }
1445 }
1446 }
1447 }
1448
1449
1450} // namespace oomph
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
char t
Definition: cfortran.h:568
Vector< double > Vector_to_symmetry_line
Vector to some point on the symmetry line along which the end of the beam is sliding.
Vector< double > Normal_to_symmetry_line
Normal vector to the symmetry line along which the end of the beam is sliding.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Fill in the element's contribution to its residual vector.
ClampedSlidingHermiteBeamBoundaryConditionElement()
Broken empty constructor.
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
Definition: shape.h:278
void initialise(const T &val)
Initialize all values in the matrix to val.
Definition: matrices.h:514
void dposition_dlagrangian_at_local_coordinate(const Vector< double > &s, DenseMatrix< double > &drdxi) const
Derivative of position vector w.r.t. the SolidFiniteElement's Lagrangian coordinates; evaluated at cu...
void get_dof_numbers_for_unknowns(std::list< std::pair< unsigned long, unsigned > > &dof_lookup_list) const
Create a list of pairs for all unknowns in this element, so that the first entry in each pair contain...
void locate_zeta(const Vector< double > &zeta, GeomObject *&geom_object_pt, Vector< double > &s, const bool &use_coordinate_as_initial_guess=false)
Find the local coordinate s in this element that corresponds to the global "intrinsic" coordinate (h...
FaceElements are elements that coincide with the faces of higher-dimensional "bulk" elements....
Definition: elements.h:4338
int & face_index()
Index of the face (a number that uniquely identifies the face in the element)
Definition: elements.h:4626
FiniteElement *& bulk_element_pt()
Pointer to higher-dimensional "bulk" element.
Definition: elements.h:4735
Vector< unsigned > Nbulk_value
A vector that will hold the number of data values at the nodes that are associated with the "bulk" el...
Definition: elements.h:4413
void get_local_coordinate_in_bulk(const Vector< double > &s, Vector< double > &s_bulk) const
Calculate the vector of local coordinate in the bulk element given the local coordinates in this Face...
Definition: elements.cc:6384
void resize_nodes(Vector< unsigned > &nadditional_data_values)
Provide additional storage for a specified number of values at the nodes of the FaceElement....
Definition: elements.h:4882
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: elements.h:4998
A general Finite Element class.
Definition: elements.h:1313
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1963
void set_nnodal_position_type(const unsigned &nposition_type)
Set the number of types required to interpolate the coordinate.
Definition: elements.h:1396
unsigned nnodal_position_type() const
Return the number of coordinate types that the element requires to interpolate the geometry between t...
Definition: elements.h:2463
void interpolated_zeta(const Vector< double > &s, Vector< double > &zeta) const
Calculate the interpolated value of zeta, the intrinsic coordinate of the element when viewed as a co...
Definition: elements.cc:4675
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
Definition: elements.cc:3962
int nodal_local_eqn(const unsigned &n, const unsigned &i) const
Return the local equation number corresponding to the i-th value at the n-th local node.
Definition: elements.h:1432
double nodal_position_gen(const unsigned &n, const unsigned &k, const unsigned &i) const
Return the value of the k-th type of the i-th positional variable at the local node n.
Definition: elements.h:2349
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
virtual void dshape_local(const Vector< double > &s, Shape &psi, DShape &dpsids) const
Function to compute the geometric shape functions and derivatives w.r.t. local coordinates at local c...
Definition: elements.h:1981
virtual void build_face_element(const int &face_index, FaceElement *face_element_pt)
Function for building a lower dimensional FaceElement on the specified face of the FiniteElement....
Definition: elements.cc:5132
double raw_dnodal_position_gen_dt(const unsigned &n, const unsigned &k, const unsigned &i) const
i-th component of time derivative (velocity) of the generalised position, dx(k,i)/dt at local node n....
Definition: elements.h:2294
unsigned nodal_dimension() const
Return the required Eulerian dimension of the nodes in this element.
Definition: elements.h:2484
double raw_nodal_position_gen(const unsigned &n, const unsigned &k, const unsigned &i) const
Return the value of the k-th type of the i-th positional variable at the local node n....
Definition: elements.h:2272
void fill_in_jacobian_from_external_by_fd(Vector< double > &residuals, DenseMatrix< double > &jacobian, const bool &fd_all_data=false)
Calculate the contributions to the jacobian from the external degrees of freedom using finite differe...
Definition: elements.cc:1199
unsigned ndof() const
Return the number of equations/dofs in the element.
Definition: elements.h:835
unsigned long eqn_number(const unsigned &ieqn_local) const
Return the global equation number corresponding to the ieqn_local-th local equation number.
Definition: elements.h:704
virtual void get_residuals(Vector< double > &residuals)
Calculate the vector of residuals of the equations in the element. By default initialise the vector t...
Definition: elements.h:980
/////////////////////////////////////////////////////////////////////
Definition: geom_objects.h:101
unsigned ndim() const
Access function to # of Eulerian coordinates.
Definition: geom_objects.h:177
virtual void d2position(const Vector< double > &zeta, RankThreeTensor< double > &ddrdzeta) const
2nd derivative of position Vector w.r.t. to coordinates: = ddrdzeta(alpha,beta,i)....
Definition: geom_objects.h:341
unsigned nlagrangian() const
Access function to # of Lagrangian coordinates.
Definition: geom_objects.h:171
Hermite Kirchhoff Love beam. Implements KirchhoffLoveBeamEquations using 2-node Hermite elements as t...
void output(std::ostream &outfile)
Output function.
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
static double Default_sigma0_value
Static default value for 2nd Piola Kirchhoff prestress.
Definition: beam_elements.h:97
static void Zero_traction_fct(const Vector< double > &xi, const Vector< double > &x, const Vector< double > &N, Vector< double > &load)
Default load function (zero traction)
GeomObject * Undeformed_beam_pt
Pointer to the GeomObject that specifies the beam's undeformed midplane.
const double & h() const
Return the non-dimensional wall thickness.
static void Unit_profile_fct(const Vector< double > &xi, const Vector< double > &x, double &h_ratio)
Default profile function (constant thickness 'h_0')
void fill_in_contribution_to_residuals_beam(Vector< double > &residuals)
Return the residuals for the equations of Kirchhoff-Love beam theory with linear constitutive equatio...
void get_non_unit_tangent(const Vector< double > &s, Vector< double > &r, Vector< double > &drds)
Get position vector to and non-unit tangent vector on wall: dr/ds.
const double & lambda_sq() const
Return the timescale ratio (non-dimensional density)
static double Default_h_value
Static default value for non-dim wall thickness.
void get_normal(const Vector< double > &s, Vector< double > &N)
Get normal vector on wall.
static double Default_lambda_sq_value
Static default value for timescale ratio (1.0 – for natural scaling)
void wall_profile(const Vector< double > &xi, const Vector< double > &x, double &h_ratio)
Get the wall profile: Pass Lagrangian & Eulerian coordinate and return the wall profile (not all of t...
const double & sigma0() const
Return the axial prestress.
void get_energy(double &pot_en, double &kin_en)
Get potential (strain) and kinetic energy of the element.
virtual void load_vector(const unsigned &intpt, const Vector< double > &xi, const Vector< double > &x, const Vector< double > &N, Vector< double > &load)
Get the load vector: Pass number of integration point (dummy), Lagr. and Eulerian coordinate and norm...
virtual void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Get FE jacobian and residuals (Jacobian done by finite differences)
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
Definition: nodes.h:1054
double value(const unsigned &i) const
Return i-th value (dofs or pinned) at this node either directly or via hanging node representation....
Definition: nodes.cc:2408
An OomphLibError object which should be thrown when an run-time error is encountered....
An OomphLibWarning object which should be created as a temporary object to issue a warning....
void shape(const Vector< double > &s, Shape &psi) const
Function to calculate the geometric shape functions at local coordinate s.
////////////////////////////////////////////////////////////////// //////////////////////////////////...
Definition: matrices.h:1370
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
double raw_lagrangian_position_gen(const unsigned &n, const unsigned &k, const unsigned &i) const
Return Generalised Lagrangian coordinate at local node n. ‘Direction’ i, ‘Type’ k....
Definition: elements.h:3897
double dshape_lagrangian(const Vector< double > &s, Shape &psi, DShape &dpsidxi) const
Calculate shape functions and derivatives w.r.t. Lagrangian coordinates at local coordinate s....
Definition: elements.cc:6710
double raw_lagrangian_position(const unsigned &n, const unsigned &i) const
Return i-th Lagrangian coordinate at local node n without using the hanging representation.
Definition: elements.h:3890
bool Solve_for_consistent_newmark_accel_flag
Flag to indicate which system of equations to solve when assigning initial conditions for time-depend...
Definition: elements.h:4302
SolidInitialCondition * Solid_ic_pt
Pointer to object that specifies the initial condition.
Definition: elements.h:4131
void fill_in_residuals_for_solid_ic(Vector< double > &residuals)
Fill in the residuals for the setup of an initial condition. The global equations are:
Definition: elements.h:4018
virtual double d2shape_lagrangian_at_knot(const unsigned &ipt, Shape &psi, DShape &dpsidxi, DShape &d2psidxi) const
Return the geometric shape functions and also first and second derivatives w.r.t. Lagrangian coordina...
Definition: elements.cc:6832
int position_local_eqn(const unsigned &n, const unsigned &k, const unsigned &j) const
Access function that returns the local equation number that corresponds to the j-th coordinate of the...
Definition: elements.h:4137
virtual double interpolated_xi(const Vector< double > &s, const unsigned &i) const
Return i-th FE-interpolated Lagrangian coordinate xi[i] at local coordinate s.
Definition: elements.cc:7104
virtual void fill_in_jacobian_from_solid_position_by_fd(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Use finite differences to calculate the Jacobian entries corresponding to the solid positions....
Definition: elements.cc:6985
void fill_in_jacobian_for_newmark_accel(DenseMatrix< double > &jacobian)
Fill in the contributions of the Jacobian matrix for the consistent assignment of the initial "accele...
Definition: elements.cc:7227
//////////////////////////////////////////////////////////////////// ////////////////////////////////...