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