shell_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: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
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 shell elements
31 #include "shell_elements.h"
32 
33 
34 namespace oomph
35 {
36 
37 //====================================================================
38 /// Default value for the Poisson ratio: 0.49 for no particular reason
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 //===================================================================
51 
52 //=======================================================================
53 /// Default load function (zero traction)
54 //=======================================================================
56  const Vector<double> &x,
57  const Vector<double>& N,
58  Vector<double>& load)
59 {
60  unsigned n_dim=load.size();
61  for (unsigned i=0;i<n_dim;i++) {load[i]=0.0;}
62 }
63 
64 
65 //=======================================================================
66 /// Static default for prestress (set to zero)
67 //=======================================================================
69 
70 
71 //======================================================================
72 /// Get normal to the shell
73 //=====================================================================
76 {
77 
78  //Set the dimension of the coordinates
79  const unsigned n_dim = 3;
80 
81  //Set the number of lagrangian coordinates
82  const unsigned n_lagrangian = 2;
83 
84  //Find out how many nodes there are
85  const unsigned n_node = nnode();
86 
87  //Find out how many positional dofs there are
88  const unsigned n_position_type = nnodal_position_type();
89 
90  //Set up memory for the shape functions
91  Shape psi(n_node,n_position_type);
92  DShape dpsidxi(n_node,n_position_type,n_lagrangian);
93 
94  //Could/Should we make this more general?
95  DShape d2psidxi(n_node,n_position_type,3);
96 
97 
98  //Call the derivatives of the shape functions:
99  // d2psidxi(i,0) = \f$ \partial^2 \psi_j / \partial^2 \xi_0^2 \f$
100  // d2psidxi(i,1) = \f$ \partial^2 \psi_j / \partial^2 \xi_1^2 \f$
101  // d2psidxi(i,2) = \f$ \partial^2 \psi_j/\partial \xi_0 \partial \xi_1 \f$
102  (void) dshape_lagrangian(s,psi,dpsidxi);
103 
104  //Calculate local values of lagrangian position and
105  //the derivative of global position wrt lagrangian coordinates
106  DenseMatrix<double> interpolated_A(2,3);
107 
108  //Initialise to zero
109  for(unsigned i=0;i<n_dim;i++)
110  {
111  for(unsigned j=0;j<n_lagrangian;j++) {interpolated_A(j,i) = 0.0;}
112  }
113 
114  //Calculate displacements and derivatives
115  for(unsigned l=0;l<n_node;l++)
116  {
117  //Loop over positional dofs
118  for(unsigned k=0;k<n_position_type;k++)
119  {
120  //Loop over displacement components (deformed position)
121  for(unsigned i=0;i<n_dim;i++)
122  {
123  double x_value = raw_nodal_position_gen(l,k,i);
124 
125  //Loop over derivative directions
126  for(unsigned j=0;j<n_lagrangian;j++)
127  {
128  interpolated_A(j,i) += x_value*dpsidxi(l,k,j);
129  }
130  }
131  }
132  }
133 
134  // Unscaled normal
135  N[0] = (interpolated_A(0,1)*interpolated_A(1,2) -
136  interpolated_A(0,2)*interpolated_A(1,1));
137  N[1] = (interpolated_A(0,2)*interpolated_A(1,0) -
138  interpolated_A(0,0)*interpolated_A(1,2));
139  N[2] = (interpolated_A(0,0)*interpolated_A(1,1) -
140  interpolated_A(0,1)*interpolated_A(1,0));
141 
142  // Normalise
143  double norm=1.0/sqrt(N[0]*N[0]+N[1]*N[1]+N[2]*N[2]);
144  for (unsigned i=0;i<3;i++) {N[i]*=norm;}
145 
146 }
147 
148 
149 //======================================================================
150 /// Get strain and bending tensors
151 //=====================================================================
153  const Vector<double>& s, DenseDoubleMatrix& strain, DenseDoubleMatrix& bend)
154  {
155  //Set the dimension of the coordinates
156  const unsigned n_dim = 3;
157 
158  //Set the number of lagrangian coordinates
159  const unsigned n_lagrangian = 2;
160 
161  //Find out how many nodes there are
162  const unsigned n_node = nnode();
163 
164  //Find out how many positional dofs there are
165  const unsigned n_position_type = nnodal_position_type();
166 
167  //Set up memory for the shape functions
168  Shape psi(n_node,n_position_type);
169  DShape dpsidxi(n_node,n_position_type,n_lagrangian);
170 
171  //Could/Should we make this more general?
172  DShape d2psidxi(n_node,n_position_type,3);
173 
174 
175  //Call the derivatives of the shape functions:
176  // d2psidxi(i,0) = \f$ \partial^2 \psi_j / \partial^2 \xi_0^2 \f$
177  // d2psidxi(i,1) = \f$ \partial^2 \psi_j / \partial^2 \xi_1^2 \f$
178  // d2psidxi(i,2) = \f$ \partial^2 \psi_j/\partial \xi_0 \partial \xi_1 \f$
179  (void) d2shape_lagrangian(s,psi,dpsidxi,d2psidxi);
180 
181  //Calculate local values of lagrangian position and
182  //the derivative of global position wrt lagrangian coordinates
184  Vector<double> accel(3,0.0);
185  DenseMatrix<double> interpolated_A(2,3);
186  DenseMatrix<double> interpolated_dAdxi(3);
187 
188  //Initialise to zero
189  for(unsigned i=0;i<n_dim;i++)
190  {
191  for(unsigned j=0;j<n_lagrangian;j++) {interpolated_A(j,i) = 0.0;}
192  for(unsigned j=0;j<3;j++) {interpolated_dAdxi(i,j) = 0.0;}
193  }
194 
195  //Calculate displacements and derivatives
196  for(unsigned l=0;l<n_node;l++)
197  {
198  //Loop over positional dofs
199  for(unsigned k=0;k<n_position_type;k++)
200  {
201  //Loop over lagrangian coordinate directions
202  for(unsigned i=0;i<n_lagrangian;i++)
203  {interpolated_xi[i] += raw_lagrangian_position_gen(l,k,i)*psi(l,k);}
204 
205  //Loop over displacement components (deformed position)
206  for(unsigned i=0;i<n_dim;i++)
207  {
208  double x_value = raw_nodal_position_gen(l,k,i);
209 
210  //Calculate interpolated (deformed) position
211  interpolated_x[i] += x_value*psi(l,k);
212 
213  //Calculate the acceleration
214  accel[i] += raw_dnodal_position_gen_dt(2,l,k,i)*psi(l,k);
215 
216  //Loop over derivative directions
217  for(unsigned j=0;j<n_lagrangian;j++)
218  {
219  interpolated_A(j,i) += x_value*dpsidxi(l,k,j);
220  }
221  //Loop over the second derivative directions
222  for(unsigned j=0;j<3;j++)
223  {
224  interpolated_dAdxi(j,i) += x_value*d2psidxi(l,k,j);
225  }
226 
227  }
228  }
229  }
230 
231  //Get the values of the undeformed parameters
232  Vector<double> interpolated_R(3);
233  DenseMatrix<double> interpolated_a(2,3);
234  RankThreeTensor<double> dadxi(2,2,3);
235 
236  //Read out the undeformed position from the geometric object
238  interpolated_a,dadxi);
239 
240  //Copy the values of the second derivatives into a slightly
241  //different data structure, where the mixed derivatives [0][1] and [1][0]
242  //are given by the single index [2]
243  DenseMatrix<double> interpolated_dadxi(3);
244  for(unsigned i=0;i<3;i++)
245  {
246  interpolated_dadxi(0,i) = dadxi(0,0,i);
247  interpolated_dadxi(1,i) = dadxi(1,1,i);
248  interpolated_dadxi(2,i) = dadxi(0,1,i);
249  }
250 
251  //Declare and calculate the undeformed and deformed metric tensor
252  double a[2][2], A[2][2], aup[2][2], Aup[2][2];
253 
254  //Assign values of A and gamma
255  for(unsigned al=0;al<2;al++)
256  {
257  for(unsigned be=0;be<2;be++)
258  {
259  //Initialise a(al,be) and A(al,be) to zero
260  a[al][be] = 0.0;
261  A[al][be] = 0.0;
262  //Now calculate the dot product
263  for(unsigned k=0;k<3;k++)
264  {
265  a[al][be] += interpolated_a(al,k)*interpolated_a(be,k);
266  A[al][be] += interpolated_A(al,k)*interpolated_A(be,k);
267  }
268  //Calculate strain tensor
269  strain(al,be) = 0.5*(A[al][be] - a[al][be]);
270  }
271  }
272 
273  //Calculate the contravariant metric tensor
274  double adet = calculate_contravariant(a,aup);
275  double Adet = calculate_contravariant(A,Aup);
276 
277  //Square roots are expensive, so let's do them once here
278  double sqrt_adet = sqrt(adet);
279  double sqrt_Adet = sqrt(Adet);
280 
281  //Calculate the normal Vectors
282  Vector<double> n(3), N(3);
283  n[0] = (interpolated_a(0,1)*interpolated_a(1,2) -
284  interpolated_a(0,2)*interpolated_a(1,1))/sqrt_adet;
285  n[1] = (interpolated_a(0,2)*interpolated_a(1,0) -
286  interpolated_a(0,0)*interpolated_a(1,2))/sqrt_adet;
287  n[2] = (interpolated_a(0,0)*interpolated_a(1,1) -
288  interpolated_a(0,1)*interpolated_a(1,0))/sqrt_adet;
289  N[0] = (interpolated_A(0,1)*interpolated_A(1,2) -
290  interpolated_A(0,2)*interpolated_A(1,1))/sqrt_Adet;
291  N[1] = (interpolated_A(0,2)*interpolated_A(1,0) -
292  interpolated_A(0,0)*interpolated_A(1,2))/sqrt_Adet;
293  N[2] = (interpolated_A(0,0)*interpolated_A(1,1) -
294  interpolated_A(0,1)*interpolated_A(1,0))/sqrt_Adet;
295 
296 
297  //Calculate the curvature tensors
298  double b[2][2], B[2][2];
299 
300  b[0][0] = n[0]*interpolated_dadxi(0,0)
301  + n[1]*interpolated_dadxi(0,1)
302  + n[2]*interpolated_dadxi(0,2);
303 
304  //Off-diagonal terms are the same
305  b[0][1] = b[1][0] = n[0]*interpolated_dadxi(2,0)
306  + n[1]*interpolated_dadxi(2,1)
307  + n[2]*interpolated_dadxi(2,2);
308 
309  b[1][1] = n[0]*interpolated_dadxi(1,0)
310  + n[1]*interpolated_dadxi(1,1)
311  + n[2]*interpolated_dadxi(1,2);
312 
313  //Deformed curvature tensor
314  B[0][0] = N[0]*interpolated_dAdxi(0,0)
315  + N[1]*interpolated_dAdxi(0,1)
316  + N[2]*interpolated_dAdxi(0,2);
317 
318  //Off-diagonal terms are the same
319  B[0][1] = B[1][0] = N[0]*interpolated_dAdxi(2,0)
320  + N[1]*interpolated_dAdxi(2,1)
321  + N[2]*interpolated_dAdxi(2,2);
322 
323  B[1][1] = N[0]*interpolated_dAdxi(1,0)
324  + N[1]*interpolated_dAdxi(1,1)
325  + N[2]*interpolated_dAdxi(1,2);
326 
327  //Set up the change of curvature tensor
328  for(unsigned i=0;i<2;i++)
329  {
330  for(unsigned j=0;j<2;j++)
331  {
332  bend(i,j) = b[i][j] - B[i][j];
333  }
334  }
335 
336  // Return undeformed and deformed determinant of in-plane metric
337  // tensor
338  return std::make_pair(adet,Adet);
339 
340  }
341 
342 
343 
344 //====================================================================
345 /// Return the residuals for the equations of KL shell theory
346 //====================================================================
349 {
350  //Set the dimension of the coordinates
351  const unsigned n_dim = 3;
352 
353  //Set the number of lagrangian coordinates
354  const unsigned n_lagrangian = 2;
355 
356  //Find out how many nodes there are
357  const unsigned n_node = nnode();
358 
359  //Find out how many positional dofs there are
360  const unsigned n_position_type = nnodal_position_type();
361 
362  //Set up memory for the shape functions
363  Shape psi(n_node,n_position_type);
364  DShape dpsidxi(n_node,n_position_type,n_lagrangian);
365 
366  //Could/Should we make this more general?
367  DShape d2psidxi(n_node,n_position_type,3);
368 
369  //Set the value of n_intpt
370  const unsigned n_intpt = integral_pt()->nweight();
371 
372  //Get Physical Variables from Element
373  //External pressure, Poisson ratio , and non-dim thickness h
374  const double nu_cached = nu();
375  const double h_cached = h();
376  const double lambda_sq_cached = lambda_sq();
377 
378  //Integers used to store the local equation numbers
379  int local_eqn=0;
380 
381  const double bending_scale = (1.0/12.0)*h_cached*h_cached;
382 
383  //Loop over the integration points
384  for(unsigned ipt=0;ipt<n_intpt;ipt++)
385  {
386  //Get the integral weight
387  double w = integral_pt()->weight(ipt);
388 
389  //Call the derivatives of the shape functions:
390  // d2psidxi(i,0) = \f$ \partial^2 \psi_j / \partial^2 \xi_0^2 \f$
391  // d2psidxi(i,1) = \f$ \partial^2 \psi_j / \partial^2 \xi_1^2 \f$
392  // d2psidxi(i,2) = \f$ \partial^2 \psi_j/\partial \xi_0 \partial \xi_1 \f$
393  double J = d2shape_lagrangian_at_knot(ipt,psi,dpsidxi,d2psidxi);
394 
395  //Premultiply the weights and the Jacobian
396  double W = w*J;
397 
398  //Calculate local values of lagrangian position and
399  //the derivative of global position wrt lagrangian coordinates
401  Vector<double> accel(3,0.0);
402  DenseMatrix<double> interpolated_A(2,3);
403  DenseMatrix<double> interpolated_dAdxi(3);
404 
405  //Initialise to zero
406  for(unsigned i=0;i<n_dim;i++)
407  {
408  for(unsigned j=0;j<n_lagrangian;j++) {interpolated_A(j,i) = 0.0;}
409  for(unsigned j=0;j<3;j++) {interpolated_dAdxi(i,j) = 0.0;}
410  }
411 
412  //Calculate displacements and derivatives
413  for(unsigned l=0;l<n_node;l++)
414  {
415  //Loop over positional dofs
416  for(unsigned k=0;k<n_position_type;k++)
417  {
418  //Loop over lagrangian coordinate directions
419  for(unsigned i=0;i<n_lagrangian;i++)
420  {interpolated_xi[i] += raw_lagrangian_position_gen(l,k,i)*psi(l,k);}
421 
422  //Loop over displacement components (deformed position)
423  for(unsigned i=0;i<n_dim;i++)
424  {
425  double x_value = raw_nodal_position_gen(l,k,i);
426 
427  //Calculate interpolated (deformed) position
428  interpolated_x[i] += x_value*psi(l,k);
429 
430  //Calculate the acceleration
431  accel[i] += raw_dnodal_position_gen_dt(2,l,k,i)*psi(l,k);
432 
433  //Loop over derivative directions
434  for(unsigned j=0;j<n_lagrangian;j++)
435  {
436  interpolated_A(j,i) += x_value*dpsidxi(l,k,j);
437  }
438  //Loop over the second derivative directions
439  for(unsigned j=0;j<3;j++)
440  {
441  interpolated_dAdxi(j,i) += x_value*d2psidxi(l,k,j);
442  }
443 
444  }
445  }
446  }
447 
448  //Get the values of the undeformed parameters
449  Vector<double> interpolated_R(3);
450  DenseMatrix<double> interpolated_a(2,3);
451  RankThreeTensor<double> dadxi(2,2,3);
452 
453  //Read out the undeformed position from the geometric object
455  interpolated_a,dadxi);
456 
457  //Copy the values of the second derivatives into a slightly
458  //different data structure, where the mixed derivatives [0][1] and [1][0]
459  //are given by the single index [2]
460  DenseMatrix<double> interpolated_dadxi(3);
461  for(unsigned i=0;i<3;i++)
462  {
463  interpolated_dadxi(0,i) = dadxi(0,0,i);
464  interpolated_dadxi(1,i) = dadxi(1,1,i);
465  interpolated_dadxi(2,i) = dadxi(0,1,i);
466  }
467 
468  //Declare and calculate the undeformed and deformed metric tensor,
469  //the strain tensor
470  double a[2][2], A[2][2], aup[2][2], Aup[2][2], gamma[2][2];
471 
472  //Assign values of A and gamma
473  for(unsigned al=0;al<2;al++)
474  {
475  for(unsigned be=0;be<2;be++)
476  {
477  //Initialise a(al,be) and A(al,be) to zero
478  a[al][be] = 0.0;
479  A[al][be] = 0.0;
480  //Now calculate the dot product
481  for(unsigned k=0;k<3;k++)
482  {
483  a[al][be] += interpolated_a(al,k)*interpolated_a(be,k);
484  A[al][be] += interpolated_A(al,k)*interpolated_A(be,k);
485  }
486  //Calculate strain tensor
487  gamma[al][be] = 0.5*(A[al][be] - a[al][be]);
488  }
489  }
490 
491  //Calculate the contravariant metric tensor
492  double adet = calculate_contravariant(a,aup);
493  double Adet = calculate_contravariant(A,Aup);
494 
495  //Square roots are expensive, so let's do them once here
496  double sqrt_adet = sqrt(adet);
497  double sqrt_Adet = sqrt(Adet);
498 
499  //Calculate the normal Vectors
500  Vector<double> n(3), N(3);
501  n[0] = (interpolated_a(0,1)*interpolated_a(1,2) -
502  interpolated_a(0,2)*interpolated_a(1,1))/sqrt_adet;
503  n[1] = (interpolated_a(0,2)*interpolated_a(1,0) -
504  interpolated_a(0,0)*interpolated_a(1,2))/sqrt_adet;
505  n[2] = (interpolated_a(0,0)*interpolated_a(1,1) -
506  interpolated_a(0,1)*interpolated_a(1,0))/sqrt_adet;
507  N[0] = (interpolated_A(0,1)*interpolated_A(1,2) -
508  interpolated_A(0,2)*interpolated_A(1,1))/sqrt_Adet;
509  N[1] = (interpolated_A(0,2)*interpolated_A(1,0) -
510  interpolated_A(0,0)*interpolated_A(1,2))/sqrt_Adet;
511  N[2] = (interpolated_A(0,0)*interpolated_A(1,1) -
512  interpolated_A(0,1)*interpolated_A(1,0))/sqrt_Adet;
513 
514 
515  //Calculate the curvature tensors
516  double b[2][2], B[2][2];
517 
518  b[0][0] = n[0]*interpolated_dadxi(0,0)
519  + n[1]*interpolated_dadxi(0,1)
520  + n[2]*interpolated_dadxi(0,2);
521 
522  //Off-diagonal terms are the same
523  b[0][1] = b[1][0] = n[0]*interpolated_dadxi(2,0)
524  + n[1]*interpolated_dadxi(2,1)
525  + n[2]*interpolated_dadxi(2,2);
526 
527  b[1][1] = n[0]*interpolated_dadxi(1,0)
528  + n[1]*interpolated_dadxi(1,1)
529  + n[2]*interpolated_dadxi(1,2);
530 
531  //Deformed curvature tensor
532  B[0][0] = N[0]*interpolated_dAdxi(0,0)
533  + N[1]*interpolated_dAdxi(0,1)
534  + N[2]*interpolated_dAdxi(0,2);
535 
536  //Off-diagonal terms are the same
537  B[0][1] = B[1][0] = N[0]*interpolated_dAdxi(2,0)
538  + N[1]*interpolated_dAdxi(2,1)
539  + N[2]*interpolated_dAdxi(2,2);
540 
541  B[1][1] = N[0]*interpolated_dAdxi(1,0)
542  + N[1]*interpolated_dAdxi(1,1)
543  + N[2]*interpolated_dAdxi(1,2);
544 
545  //Set up the change of curvature tensor
546  double kappa[2][2];
547 
548  for(unsigned i=0;i<2;i++)
549  {
550  for(unsigned j=0;j<2;j++)
551  {
552  kappa[i][j] = b[i][j] - B[i][j];
553  }
554  }
555 
556  //Calculate the plane stress stiffness tensor
557  double Et[2][2][2][2];
558 
559  //Premultiply some constants
560  //NOTE C2 has 1.0-Nu in the denominator
561  //double C1 = 1.0/(2.0*(1.0+nu_cached)), C2 = 2.0*nu_cached/(1.0-nu_cached);
562 
563  // Some constants
564  double C1=0.5*(1.0-nu_cached);
565  double C2=nu_cached;
566 
567  //Loop over first index
568  for(unsigned i=0;i<2;i++)
569  {
570  for(unsigned j=0;j<2;j++)
571  {
572  for(unsigned k=0;k<2;k++)
573  {
574  for(unsigned l=0;l<2;l++)
575  {
576  //Et[i][j][k][l] = C1*(aup[i][k]*aup[j][l] + aup[i][l]*aup[j][k]
577  // + C2*aup[i][j]*aup[k][l]);
578  Et[i][j][k][l] = C1*(aup[i][k]*aup[j][l] + aup[i][l]*aup[j][k])
579  + C2*aup[i][j]*aup[k][l];
580  }
581  }
582  }
583  }
584 
585  //Define the load Vector
586  Vector<double> f(3);
587  //Determine the load
588  load_vector(ipt,interpolated_xi,interpolated_x,N,f);
589 
590  //Scale the load by the bending stiffness scale
591  //for(unsigned i=0;i<3;i++){f[i] *= bending_scale/(1.0-nu_cached*nu_cached);}
592 
593  //Allocate storage for variations of the the normal vector
594  double normal_var[3][3];
595 
596 //=====EQUATIONS OF ELASTICITY FROM PRINCIPLE OF VIRTUAL DISPLACEMENTS========
597 
598  //Little tensor to handle the mixed derivative terms
599  unsigned mix[2][2] = {{0,2},{2,1}};
600 
601  //Loop over the number of nodes
602  for(unsigned l=0;l<n_node;l++)
603  {
604  //Loop over the type of degree of freedom
605  for(unsigned k=0;k<n_position_type;k++)
606  {
607  //Setup components of variations in the normal vector
608  normal_var[0][0] = 0.0;
609  normal_var[0][1] = (-interpolated_A(1,2)*dpsidxi(l,k,0)
610  + interpolated_A(0,2)*dpsidxi(l,k,1))/sqrt_Adet;
611  normal_var[0][2] = (interpolated_A(1,1)*dpsidxi(l,k,0)
612  - interpolated_A(0,1)*dpsidxi(l,k,1))/sqrt_Adet;
613 
614  normal_var[1][0] = (interpolated_A(1,2)*dpsidxi(l,k,0)
615  - interpolated_A(0,2)*dpsidxi(l,k,1))/sqrt_Adet;
616  normal_var[1][1] = 0.0;
617  normal_var[1][2] = (-interpolated_A(1,0)*dpsidxi(l,k,0)
618  + interpolated_A(0,0)*dpsidxi(l,k,1))/sqrt_Adet;
619 
620  normal_var[2][0] = (-interpolated_A(1,1)*dpsidxi(l,k,0)
621  + interpolated_A(0,1)*dpsidxi(l,k,1))/sqrt_Adet;
622  normal_var[2][1] = (interpolated_A(1,0)*dpsidxi(l,k,0)
623  - interpolated_A(0,0)*dpsidxi(l,k,1))/sqrt_Adet;
624  normal_var[2][2] = 0.0;
625 
626  //Loop over the coordinate direction
627  for(unsigned i=0;i<3;i++)
628  {
629  //Get the local equation number
630  local_eqn = position_local_eqn(l,k,i);
631 
632  //If it's not a boundary condition
633  if(local_eqn >= 0)
634  {
635  //Temporary to hold a common part of the bending terms
636  double bending_dot_product_multiplier =
637  ((A[1][1]*interpolated_A(0,i) -
638  A[1][0]*interpolated_A(1,i))*dpsidxi(l,k,0)
639  + (A[0][0]*interpolated_A(1,i) -
640  A[0][1]*interpolated_A(0,i))*dpsidxi(l,k,1))/Adet;
641 
642  //Combine the cross and dot product bending terms
643  for(unsigned i2=0;i2<3;i2++)
644  {normal_var[i][i2] -= N[i2]*bending_dot_product_multiplier;}
645 
646  //Add in external forcing
647  residuals[local_eqn] -= f[i]/h_cached*psi(l,k)*W*sqrt_Adet;
648 
649  //Storage for all residual terms that are multipled by the
650  //Jacobian and area of the undeformed shell
651  //Start with the acceleration term
652  double other_residual_terms = lambda_sq_cached*accel[i]*psi(l,k);
653 
654  //Now need the loops over the Greek indicies
655  for(unsigned al=0;al<2;al++)
656  {
657  for(unsigned be=0;be<2;be++)
658  {
659  // Membrane prestress
660  other_residual_terms +=
661  *Prestress_pt(al,be)
662  *interpolated_A(al,i)
663  *dpsidxi(l,k,be);
664 
665  // Remaining terms
666  for(unsigned ga=0;ga<2;ga++)
667  {
668  for(unsigned de=0;de<2;de++)
669  {
671  {
672  //Pure membrane term
673  other_residual_terms +=
674  Et[al][be][ga][de]*gamma[al][be]
675  *interpolated_A(ga,i)
676  *dpsidxi(l,k,de);
677  }
678 
679  //Bending terms
680  other_residual_terms -=
681  bending_scale*Et[al][be][ga][de]*kappa[al][be]*
682  (N[i]*d2psidxi(l,k,mix[ga][de])
683  //Cross and dot product terms
684  + normal_var[i][0]*interpolated_dAdxi(mix[ga][de],0)
685  + normal_var[i][1]*interpolated_dAdxi(mix[ga][de],1)
686  + normal_var[i][2]*interpolated_dAdxi(mix[ga][de],2));
687  }
688  }
689  }
690  }
691 
692  residuals[local_eqn] += other_residual_terms*W*sqrt_adet;
693  }
694  }
695  }
696  }
697 
698  } //End of loop over the integration points
699 
700 }
701 
702 //=========================================================================
703 /// Return the jacobian is calculated by finite differences by default,
704 //=========================================================================
707  DenseMatrix<double> &jacobian)
708 {
709  //Call the element's residuals vector
711 
712  //Solve for the consistent acceleration in Newmark scheme?
714  {
716  return;
717  }
718 
719  //Allocate storage for the full residuals of the element
720  unsigned n_dof = ndof();
721  Vector<double> full_residuals(n_dof);
722 
723 //Call the full residuals
724  get_residuals(full_residuals);
725 
726  //Get the solid entries in the jacobian using finite differences
728  fill_in_jacobian_from_solid_position_by_fd(full_residuals,jacobian);
729 
730  //Get the entries from the external data, usually load terms
732  fill_in_jacobian_from_external_by_fd(full_residuals,jacobian);
733 }
734 
735 //=======================================================================
736 /// Compute the potential (strain) and kinetic energy of the
737 /// element.
738 //=======================================================================
739 void KirchhoffLoveShellEquations::get_energy(double& pot_en, double& kin_en)
740 {
741  // Initialise
742  pot_en=0.0;
743  kin_en=0.0;
744 
745  //Set the dimension of the coordinates
746  const unsigned n_dim = Undeformed_midplane_pt->ndim();
747 
748  //Set the number of lagrangian coordinates
749  const unsigned n_lagrangian = Undeformed_midplane_pt->nlagrangian();
750 
751  //Find out how many nodes there are
752  const unsigned n_node = nnode();
753 
754  //Find out how many positional dofs there are
755  const unsigned n_position_type = nnodal_position_type();
756 
757  //Set up memory for the shape functions
758  Shape psi(n_node,n_position_type);
759  DShape dpsidxi(n_node,n_position_type,n_lagrangian);
760  DShape d2psidxi(n_node,n_position_type,3);
761 
762  //Set the value of n_intpt
763  const unsigned n_intpt = integral_pt()->nweight();
764 
765  //Get Physical Variables from Element
766  //External pressure, Poisson ratio , and non-dim thickness h
767  const double nu_cached = nu();
768  const double h_cached = h();
769  const double lambda_sq_cached = lambda_sq();
770 
771 #ifdef PARANOID
772  //Check for non-zero prestress
773  if( (prestress(0,0)!=0) || (prestress(1,0)!=0) || (prestress(1,1)!=0) )
774  {
775  std::string error_message =
776  "Warning: Not sure if the energy is computed correctly\n";
777  error_message += " for nonzero prestress\n";
778 
779  oomph_info << error_message << std::endl;
780 
781 // throw OomphLibWarning(error_message,
782 // "KirchhoffLoveShellEquations::get_energy()",
783 // OOMPH_EXCEPTION_LOCATION);
784  }
785 #endif
786 
787  // Setup storage for various quantities
788  Vector<double> veloc(n_dim);
789  Vector<double> interpolated_xi(n_lagrangian);
790  DenseMatrix<double> interpolated_A(n_lagrangian,n_dim);
791  DenseMatrix<double> interpolated_dAdxi(3);
792 
793  const double bending_scale = (1.0/12.0)*h_cached*h_cached;
794 
795  //Loop over the integration points
796  for(unsigned ipt=0;ipt<n_intpt;ipt++)
797  {
798  //Get the integral weight
799  double w = integral_pt()->weight(ipt);
800 
801  //Call the derivatives of the shape functions:
802  double J = d2shape_lagrangian_at_knot(ipt,psi,dpsidxi,d2psidxi);
803 
804  //Premultiply the weights and the Jacobian
805  double W = w*J;
806 
807  //Initialise values to zero
808  for(unsigned i=0;i<n_dim;i++)
809  {
810  veloc[i]=0.0;
811  for(unsigned j=0;j<n_lagrangian;j++)
812  {
813  interpolated_A(j,i) = 0.0;
814  }
815  for(unsigned j=0;j<3;j++)
816  {
817  interpolated_dAdxi(i,j) = 0.0;
818  }
819  }
820  for(unsigned j=0;j<n_lagrangian;j++)
821  {
822  interpolated_xi[j] = 0.0;
823  }
824 
825  //Calculate velocity and derivatives
826  for(unsigned l=0;l<n_node;l++)
827  {
828  //Loop over positional dofs
829  for(unsigned k=0;k<n_position_type;k++)
830  {
831  //Loop over lagrangian coordinate directions
832  for(unsigned i=0;i<n_lagrangian;i++)
833  {
834  interpolated_xi[i] += lagrangian_position_gen(l,k,i)*psi(l,k);
835  }
836 
837  //Loop over displacement components
838  for(unsigned i=0;i<n_dim;i++)
839  {
840  veloc[i] += dnodal_position_gen_dt(1,l,k,i)*psi(l,k);
841 
842  double x_value = nodal_position_gen(l,k,i);
843 
844  //Loop over derivative directions
845  for(unsigned j=0;j<n_lagrangian;j++)
846  {
847  interpolated_A(j,i) += x_value*dpsidxi(l,k,j);
848  }
849  //Loop over the second derivative directions
850  for(unsigned j=0;j<3;j++)
851  {
852  interpolated_dAdxi(j,i) += x_value*d2psidxi(l,k,j);
853  }
854  }
855  }
856  }
857 
858  // Get square of veloc
859  double veloc_sq=0;
860  for(unsigned i=0;i<n_dim;i++)
861  {
862  veloc_sq += veloc[i]*veloc[i];
863  }
864 
865  //Get the values of the undeformed parameters
866  Vector<double> interpolated_r(n_dim);
867  DenseMatrix<double> interpolated_a(n_lagrangian,n_dim);
868  RankThreeTensor<double> dadxi(n_lagrangian,n_lagrangian,n_dim);
869 
870  //Read out the undeformed position from the geometric object
871  Undeformed_midplane_pt->d2position(interpolated_xi,
872  interpolated_r,
873  interpolated_a,
874  dadxi);
875 
876  //Copy the values of the second derivatives into a slightly
877  //different data structure, where the mixed derivatives [0][1] and [1][0]
878  //are given by the single index [2]
879  DenseMatrix<double> interpolated_dadxi(3);
880  for(unsigned i=0;i<3;i++)
881  {
882  interpolated_dadxi(0,i) = dadxi(0,0,i);
883  interpolated_dadxi(1,i) = dadxi(1,1,i);
884  interpolated_dadxi(2,i) = dadxi(0,1,i);
885  }
886 
887  //Declare and calculate the undeformed and deformed metric tensor,
888  //the strain tensor
889  double a[2][2], A[2][2], aup[2][2], Aup[2][2], gamma[2][2];
890 
891  //Assign values of A and gamma
892  for(unsigned al=0;al<2;al++)
893  {
894  for(unsigned be=0;be<2;be++)
895  {
896  //Initialise a(al,be) and A(al,be) to zero
897  a[al][be] = 0.0;
898  A[al][be] = 0.0;
899  //Now calculate the dot product
900  for(unsigned k=0;k<n_dim;k++)
901  {
902  a[al][be] += interpolated_a(al,k)*interpolated_a(be,k);
903  A[al][be] += interpolated_A(al,k)*interpolated_A(be,k);
904  }
905  //Calculate strain tensor
906  gamma[al][be] = 0.5*(A[al][be] - a[al][be]);
907  }
908  }
909 
910  //Calculate the contravariant metric tensor
911  double adet = calculate_contravariant(a,aup);
912  double Adet = calculate_contravariant(A,Aup);
913 
914  //Square roots are expensive, so let's do them once here
915  double sqrt_adet = sqrt(adet);
916  double sqrt_Adet = sqrt(Adet);
917 
918  //Calculate the normal Vectors
919  Vector<double> n(3), N(3);
920  n[0] = (interpolated_a(0,1)*interpolated_a(1,2) -
921  interpolated_a(0,2)*interpolated_a(1,1))/sqrt_adet;
922  n[1] = (interpolated_a(0,2)*interpolated_a(1,0) -
923  interpolated_a(0,0)*interpolated_a(1,2))/sqrt_adet;
924  n[2] = (interpolated_a(0,0)*interpolated_a(1,1) -
925  interpolated_a(0,1)*interpolated_a(1,0))/sqrt_adet;
926  N[0] = (interpolated_A(0,1)*interpolated_A(1,2) -
927  interpolated_A(0,2)*interpolated_A(1,1))/sqrt_Adet;
928  N[1] = (interpolated_A(0,2)*interpolated_A(1,0) -
929  interpolated_A(0,0)*interpolated_A(1,2))/sqrt_Adet;
930  N[2] = (interpolated_A(0,0)*interpolated_A(1,1) -
931  interpolated_A(0,1)*interpolated_A(1,0))/sqrt_Adet;
932 
933 
934  //Calculate the curvature tensors
935  double b[2][2], B[2][2];
936 
937  b[0][0] = n[0]*interpolated_dadxi(0,0)
938  + n[1]*interpolated_dadxi(0,1)
939  + n[2]*interpolated_dadxi(0,2);
940 
941  //Off-diagonal terms are the same
942  b[0][1] = b[1][0] = n[0]*interpolated_dadxi(2,0)
943  + n[1]*interpolated_dadxi(2,1)
944  + n[2]*interpolated_dadxi(2,2);
945 
946  b[1][1] = n[0]*interpolated_dadxi(1,0)
947  + n[1]*interpolated_dadxi(1,1)
948  + n[2]*interpolated_dadxi(1,2);
949 
950  //Deformed curvature tensor
951  B[0][0] = N[0]*interpolated_dAdxi(0,0)
952  + N[1]*interpolated_dAdxi(0,1)
953  + N[2]*interpolated_dAdxi(0,2);
954 
955  //Off-diagonal terms are the same
956  B[0][1] = B[1][0] = N[0]*interpolated_dAdxi(2,0)
957  + N[1]*interpolated_dAdxi(2,1)
958  + N[2]*interpolated_dAdxi(2,2);
959 
960  B[1][1] = N[0]*interpolated_dAdxi(1,0)
961  + N[1]*interpolated_dAdxi(1,1)
962  + N[2]*interpolated_dAdxi(1,2);
963 
964  //Set up the change of curvature tensor
965  double kappa[2][2];
966 
967  for(unsigned i=0;i<2;i++)
968  {
969  for(unsigned j=0;j<2;j++)
970  {
971  kappa[i][j] = b[i][j] - B[i][j];
972  }
973  }
974 
975  //Calculate the plane stress stiffness tensor
976  double Et[2][2][2][2];
977 
978  // Some constants
979  double C1=0.5*(1.0-nu_cached);
980  double C2=nu_cached;
981 
982  //Loop over first index
983  for(unsigned i=0;i<2;i++)
984  {
985  for(unsigned j=0;j<2;j++)
986  {
987  for(unsigned k=0;k<2;k++)
988  {
989  for(unsigned l=0;l<2;l++)
990  {
991  Et[i][j][k][l] = C1*(aup[i][k]*aup[j][l] + aup[i][l]*aup[j][k])
992  + C2*aup[i][j]*aup[k][l];
993  }
994  }
995  }
996  }
997 
998  // Add contributions to potential energy
999  double pot_en_cont = 0.0;
1000 
1001  for(unsigned i=0;i<2;i++)
1002  {
1003  for(unsigned j=0;j<2;j++)
1004  {
1005  for(unsigned k=0;k<2;k++)
1006  {
1007  for(unsigned l=0;l<2;l++)
1008  {
1009  pot_en_cont+= Et[i][j][k][l] *
1010  (gamma[i][j]*gamma[k][l] + bending_scale*kappa[i][j]*kappa[k][l]);
1011  }
1012  }
1013  }
1014  }
1015  pot_en += pot_en_cont*0.5*W*sqrt_adet;
1016 
1017  // Add contribution to kinetic energy
1018  kin_en+=0.5*lambda_sq_cached*veloc_sq*W*sqrt_adet;
1019 
1020  } //End of loop over the integration points
1021 }
1022 
1023 
1024 
1025 //===================================================================
1026 /// Get integral of instantaneous rate of work done on
1027 /// the wall due to the load returned by the virtual
1028 /// function load_vector_for_rate_of_work_computation(...).
1029 /// In the current class
1030 /// the latter function simply de-references the external
1031 /// load but this can be overloaded in derived classes
1032 /// (e.g. in FSI elements) to determine the rate of work done
1033 /// by individual constituents of this load (e.g. the fluid load
1034 /// in an FSI problem).
1035 //===================================================================
1037 {
1038  // Initialise
1039  double rate_of_work_integral=0.0;
1040 
1041  // The number of dimensions
1042  const unsigned n_dim = 3;
1043 
1044  // The number of lagrangian coordinates
1045  const unsigned n_lagrangian = 2;
1046 
1047  // The number of nodes
1048  const unsigned n_node = nnode();
1049 
1050  // Number of integration points
1051  unsigned n_intpt = integral_pt()->nweight();
1052 
1053  //Find out how many positional dofs there are
1054  unsigned n_position_type = nnodal_position_type();
1055 
1056  // Vector to hold local coordinate in element
1057  Vector<double> s(n_lagrangian);
1058 
1059  // Set up storage for shape functions and derivatives of shape functions
1060  Shape psi(n_node, n_position_type);
1061  DShape dpsidxi(n_node, n_position_type, n_lagrangian);
1062 
1063  // Storage for jacobian of mapping from local to Lagrangian coords
1064  DenseMatrix<double> jacobian(n_lagrangian);
1065 
1066  // Storage for velocity vector
1067  Vector<double> velocity(n_dim);
1068 
1069  // Storage for load vector
1070  Vector<double> load(n_dim);
1071 
1072  // Storage for normal vector
1073  Vector<double> normal(n_dim);
1074 
1075  // Storage for Lagrangian and Eulerian coordinates
1076  Vector<double> xi(n_lagrangian);
1077  Vector<double> x(n_dim);
1078 
1079  // Storage for covariant vectors
1080  DenseMatrix<double> interpolated_A(n_lagrangian,n_dim);
1081 
1082  //Loop over the integration points
1083  for (unsigned ipt=0;ipt<n_intpt;ipt++)
1084  {
1085  //Get the integral weight
1086  double w = integral_pt()->weight(ipt);
1087 
1088  // Get local coords at knot
1089  for(unsigned i=0;i<n_lagrangian;i++)
1090  {
1091  s[i] = integral_pt()->knot(ipt,i);
1092  }
1093 
1094  // Get shape functions, derivatives, determinant of Jacobian of mapping
1095  double J = dshape_lagrangian(s,psi,dpsidxi);
1096 
1097  // Premultiply the weights and the Jacobian
1098  double W = w*J;
1099 
1100  // Initialise velocity, x and interpolated_A to zero
1101  for(unsigned i=0;i<n_dim;i++)
1102  {
1103  velocity[i]=0.0;
1104  x[i]=0.0;
1105  for(unsigned j=0;j<n_lagrangian;j++)
1106  {
1107  interpolated_A(j,i) = 0.0;
1108  }
1109  }
1110 
1111  // Initialise xi to zero
1112  for(unsigned i=0;i<n_lagrangian;i++)
1113  {
1114  xi[i]=0.0;
1115  }
1116 
1117  // Calculate velocity, coordinates and derivatives
1118  for(unsigned l=0;l<n_node;l++)
1119  {
1120  // Loop over positional dofs
1121  for(unsigned k=0;k<n_position_type;k++)
1122  {
1123  // Loop over lagrangian coordinate directions
1124  for(unsigned i=0;i<n_lagrangian;i++)
1125  {
1126  xi[i] += lagrangian_position_gen(l,k,i)*psi(l,k);
1127  }
1128 
1129  // Loop over displacement components
1130  for(unsigned i=0;i<n_dim;i++)
1131  {
1132  velocity[i] += dnodal_position_gen_dt(1,l,k,i)*psi(l,k);
1133 
1134  double x_value = nodal_position_gen(l,k,i);
1135  x[i] += x_value*psi(l,k);
1136 
1137  // Loop over derivative directions
1138  for(unsigned j=0;j<n_lagrangian;j++)
1139  {
1140  interpolated_A(j,i) += x_value*dpsidxi(l,k,j);
1141  }
1142  }
1143  }
1144  }
1145 
1146  // Get deformed metric tensor
1147  double A[2][2];
1148  for(unsigned al=0;al<2;al++)
1149  {
1150  for(unsigned be=0;be<2;be++)
1151  {
1152  // Initialise A(al,be) to zero
1153  A[al][be] = 0.0;
1154  // Calculate the dot product
1155  for(unsigned k=0;k<n_dim;k++)
1156  {
1157  A[al][be] += interpolated_A(al,k)*interpolated_A(be,k);
1158  }
1159  }
1160  }
1161 
1162  // Get determinant of metric tensor
1163  double A_det = A[0][0]*A[1][1] - A[0][1]*A[1][0];
1164 
1165  // Get outer unit normal
1166  get_normal(s,normal);
1167 
1168  // Get load vector
1169  load_vector_for_rate_of_work_computation(ipt, xi, x, normal,load);
1170 
1171  // Local rate of work:
1172  double rate_of_work=0.0;
1173  for (unsigned i=0;i<n_dim;i++)
1174  {
1175  rate_of_work+=load[i]*velocity[i];
1176  }
1177 
1178  // Add rate of work
1179  rate_of_work_integral+=rate_of_work/h()*W*A_det;
1180  }
1181 
1182  return rate_of_work_integral;
1183 }
1184 
1185 
1186 
1187 
1188 //===================================================================
1189 /// The output function: position, veloc and accel.
1190 //===================================================================
1192  const unsigned &n_plot)
1193 {
1194  //Local variables
1195  Vector<double> s(2);
1196 
1197  //Tecplot header info
1198  outfile << "ZONE I=" << n_plot << ", J=" << n_plot << std::endl;
1199 
1200  //Loop over plot points
1201  for(unsigned l2=0;l2<n_plot;l2++)
1202  {
1203  s[1] = -1.0 + l2*2.0/(n_plot-1);
1204  for(unsigned l1=0;l1<n_plot;l1++)
1205  {
1206  s[0] = -1.0 + l1*2.0/(n_plot-1);
1207 
1208  //Output
1209  outfile << interpolated_x(s,0) << " "
1210  << interpolated_x(s,1) << " "
1211  << interpolated_x(s,2) << " "
1212  << interpolated_dxdt(s,0,1) << " "
1213  << interpolated_dxdt(s,1,1) << " "
1214  << interpolated_dxdt(s,2,1) << " "
1215  << interpolated_dxdt(s,0,2) << " "
1216  << interpolated_dxdt(s,1,2) << " "
1217  << interpolated_dxdt(s,2,2) << " "
1218  << std::endl;
1219  }
1220  }
1221  outfile << std::endl;
1222 }
1223 
1224 
1225 //===================================================================
1226 /// The output function
1227 //===================================================================
1228 void HermiteShellElement::output(std::ostream &outfile, const unsigned &n_plot)
1229 {
1230  //Local variables
1231  Vector<double> s(2);
1232 
1233  //Tecplot header info
1234  outfile << "ZONE I=" << n_plot << ", J=" << n_plot << std::endl;
1235 
1236  //Loop over plot points
1237  for(unsigned l2=0;l2<n_plot;l2++)
1238  {
1239  s[1] = -1.0 + l2*2.0/(n_plot-1);
1240  for(unsigned l1=0;l1<n_plot;l1++)
1241  {
1242  s[0] = -1.0 + l1*2.0/(n_plot-1);
1243 
1244  //Output the x and y positions
1245  outfile << interpolated_x(s,0) << " " << interpolated_x(s,1) << " "
1246  << interpolated_x(s,2) << " ";
1247  outfile << interpolated_xi(s,0) << " " << interpolated_xi(s,1)
1248  << std::endl;
1249  }
1250  }
1251  outfile << std::endl;
1252 }
1253 
1254 
1255 //===================================================================
1256 /// The output function
1257 //===================================================================
1258 void HermiteShellElement::output(FILE* file_pt, const unsigned &n_plot)
1259 {
1260  //Local variables
1261  Vector<double> s(2);
1262 
1263  //Tecplot header info
1264  fprintf(file_pt,"ZONE I=%i, J=%i \n",n_plot,n_plot);
1265 
1266  //Loop over plot points
1267  for(unsigned l2=0;l2<n_plot;l2++)
1268  {
1269  s[1] = -1.0 + l2*2.0/(n_plot-1);
1270  for(unsigned l1=0;l1<n_plot;l1++)
1271  {
1272  s[0] = -1.0 + l1*2.0/(n_plot-1);
1273 
1274  //Output the x and y positions
1275  fprintf(file_pt,"%g %g %g ",
1276  interpolated_x(s,0),
1277  interpolated_x(s,1),
1278  interpolated_x(s,2));
1279  fprintf(file_pt,"%g %g \n ",
1280  interpolated_xi(s,0),
1281  interpolated_xi(s,1));
1282  }
1283  }
1284  fprintf(file_pt,"\n");
1285 }
1286 
1287 
1288 
1289 
1290 //=========================================================================
1291 /// Define the dposition function. This is used to set no-slip boundary
1292 /// conditions in some FSI problems.
1293 //=========================================================================
1296  const Vector<double> &s, DenseMatrix<double> &drdxi) const
1297 {
1298 #ifdef PARANOID
1299  if(Undeformed_midplane_pt==0)
1300  {
1301  throw
1302  OomphLibError(
1303  "Undeformed_midplane_pt has not been set",
1304  "FSIDiagHermiteShellElement::dposition_dlagrangian_at_local_coordinate()",
1305  OOMPH_EXCEPTION_LOCATION);
1306  }
1307 #endif
1308 
1309  //Find the dimension of the global coordinate
1310  unsigned n_dim = Undeformed_midplane_pt->ndim();
1311 
1312  //Find the number of lagrangian coordinates
1313  unsigned n_lagrangian = Undeformed_midplane_pt->nlagrangian();
1314 
1315  //Find out how many nodes there are
1316  unsigned n_node = nnode();
1317 
1318  //Find out how many positional dofs there are
1319  unsigned n_position_type = this->nnodal_position_type();
1320 
1321  //Set up memory for the shape functions
1322  Shape psi(n_node,n_position_type);
1323  DShape dpsidxi(n_node,n_position_type,n_lagrangian);
1324 
1325  //Get the derivatives of the shape functions w.r.t local coordinates
1326  //at this point
1327  dshape_lagrangian(s,psi,dpsidxi);
1328 
1329  //Initialise the derivatives to zero
1330  drdxi.initialise(0.0);
1331 
1332  //Loop over the nodes
1333  for(unsigned l=0;l<n_node;l++)
1334  {
1335  //Loop over the positional types
1336  for(unsigned k=0;k<n_position_type;k++)
1337  {
1338  //Loop over the dimensions
1339  for(unsigned i=0;i<n_dim;i++)
1340  {
1341  //Loop over the lagrangian coordinates
1342  for(unsigned j=0;j<n_lagrangian;j++)
1343  {
1344  //Add the contribution to the derivative
1345  drdxi(j,i) += nodal_position_gen(l,k,i)*dpsidxi(l,k,j);
1346  }
1347  }
1348  }
1349  }
1350 }
1351 
1352 
1353 //=============================================================================
1354 /// Create a list of pairs for all unknowns in this element,
1355 /// so that the first entry in each pair contains the global equation
1356 /// number of the unknown, while the second one contains the number
1357 /// of the "DOF types" that this unknown is associated with.
1358 /// (Function can obviously only be called if the equation numbering
1359 /// scheme has been set up.)
1360 /// This element is only in charge of the solid dofs.
1361 //=============================================================================
1363  std::list<std::pair<unsigned long,unsigned> >& dof_lookup_list) const
1364 {
1365 
1366  // temporary pair (used to store dof lookup prior to being added to list)
1367  std::pair<unsigned long,unsigned> dof_lookup;
1368 
1369  // number of nodes
1370  const unsigned n_node = this->nnode();
1371 
1372  //Get the number of position dofs and dimensions at the node
1373  const unsigned n_position_type = nnodal_position_type();
1374  const unsigned nodal_dim = nodal_dimension();
1375 
1376  //Integer storage for local unknown
1377  int local_unknown=0;
1378 
1379  //Loop over the nodes
1380  for(unsigned n=0;n<n_node;n++)
1381  {
1382  //Loop over position dofs
1383  for(unsigned k=0;k<n_position_type;k++)
1384  {
1385  //Loop over dimension
1386  for(unsigned i=0;i<nodal_dim;i++)
1387  {
1388  //If the variable is free
1389  local_unknown = position_local_eqn(n,k,i);
1390 
1391  // ignore pinned values
1392  if (local_unknown >= 0)
1393  {
1394  // store dof lookup in temporary pair: First entry in pair
1395  // is global equation number; second entry is dof type
1396  dof_lookup.first = this->eqn_number(local_unknown);
1397  dof_lookup.second = 0;
1398 
1399  // add to list
1400  dof_lookup_list.push_front(dof_lookup);
1401 
1402  }
1403  }
1404  }
1405  }
1406 }
1407 
1408 
1409 ////////////////////////////////////////////////////////////////////////
1410 ////////////////////////////////////////////////////////////////////////
1411 ////////////////////////////////////////////////////////////////////////
1412 
1413 
1414 //===========================================================================
1415 /// Constructor, takes the pointer to the "bulk" element, the
1416 /// index of the fixed local coordinate and its value represented
1417 /// by an integer (+/- 1), indicating that the face is located
1418 /// at the max. or min. value of the "fixed" local coordinate
1419 /// in the bulk element.
1420 //===========================================================================
1423  FiniteElement* const &bulk_el_pt,
1424  const int &face_index) :
1426 {
1427 
1428  // Let the bulk element build the FaceElement, i.e. setup the pointers
1429  // to its nodes (by referring to the appropriate nodes in the bulk
1430  // element), etc.
1431  bulk_el_pt->build_face_element(face_index,this);
1432 
1433  // Overwrite default assignments:
1434 
1435  // Number of position types (for automatic local equation numbering)
1436  // must be re-set to the shell case because the constraint
1437  // equation makes contributions to all 3 (coordinate directions) x
1438  // 4 (types of positional dofs) = 12 dofs associated with each node in the
1439  // shell element. But keep reading...
1441 
1442 
1443  // The trouble with the above is that the parametrisation of the
1444  // element geometry (used in J_eulerian and elsewhere)
1445  // is based on the 4 (2 nodes x 2 types -- value and slope) shape functions
1446  // in the 1D element. In general this is done by looping over the
1447  // nodal position types (2 per node) and referring to the
1448  // relevant position Data at the nodes indirectly via the
1449  // bulk_position_type(...) function. This is now screwed up
1450  // because we've bumped up the number of nodal position types
1451  // to four! The slightly hacky way around this is resize
1452  // the lookup scheme underlying the bulk_position_type(...)
1453  // to four types of degree of freedom per node so it's consistent
1454  // with the above assignement. However, we only have shape
1455  // functions associated with the two types of degree of freedom,
1456  // so (below) we overload shape(...) and dshape_local(...)
1457  // to that it returns zero for all "superfluous" shape functions.
1459 
1460  // Make some dummy assigments to the position types that
1461  // correspond to the two superfluous ones -- the values associated
1462  // with these will never be used since we're setting the
1463  // associated shape functions to zero!
1464  bulk_position_type(2)=0;
1465  bulk_position_type(3)=0;
1466 
1467  // Resize normal to clamping plane and initialise for clamping
1468  // being applied in a plane where z=const.
1469  Normal_to_clamping_plane.resize(3);
1470  Normal_to_clamping_plane[0]=0.0;
1471  Normal_to_clamping_plane[1]=0.0;
1472  Normal_to_clamping_plane[2]=1.0;
1473 
1474  // We need two additional values (for the two types of the Lagrange
1475  // multiplier interpolated by the 1D Hermite basis functions)
1476  // at the element's two nodes
1477  Vector<unsigned> nadditional_data_values(2);
1478  nadditional_data_values[0]=2;
1479  nadditional_data_values[1]=2;
1480  resize_nodes(nadditional_data_values);
1481 
1482 }
1483 
1484 
1485 
1486 //===========================================================================
1487 /// Add the element's contribution to its residual vector
1488 //===========================================================================
1491 {
1492  // Get position vector to and normal vector on wall (from bulk element)
1493  HermiteShellElement* bulk_el_pt=
1494  dynamic_cast<HermiteShellElement*>(bulk_element_pt());
1495 
1496  // Local coordinates in bulk
1497  Vector<double> s_bulk(2);
1498 
1499  // Local coordinate in FaceElement:
1500  Vector<double> s(1);
1501 
1502  // Normal to shell
1503  Vector<double> shell_normal(3);
1504  Vector<double> shell_normal_plus(3);
1505 
1506  //Find out how many nodes there are
1507  const unsigned n_node = nnode();
1508 
1509  // Types of (nonzero) shape/test fcts (for Lagrange multiplier) in
1510  // this element
1511  unsigned n_type=2;
1512 
1513  //Integer to store the local equation number
1514  int local_eqn=0;
1515 
1516  //Set up memory for the shape functions:
1517 
1518  // Shape fcts
1519  Shape psi(n_node,n_type);
1520 
1521  //Set # of integration points
1522  const unsigned n_intpt = integral_pt()->nweight();
1523 
1524  //Loop over the integration points
1525  for(unsigned ipt=0;ipt<n_intpt;ipt++)
1526  {
1527  //Get the integral weight
1528  double w = integral_pt()->weight(ipt);
1529 
1530  // Get integration point in local coordinates
1531  s[0]=integral_pt()->knot(ipt,0);
1532 
1533  // Get shape functions
1534  shape(s,psi);
1535 
1536  // Get jacobian
1537  double J=J_eulerian(s);
1538 
1539  // Assemble the Lagrange multiplier
1540  double lambda=0.0;
1541  for(unsigned j=0;j<n_node;j++)
1542  {
1543  for(unsigned k=0;k<n_type;k++)
1544  {
1545  // The (additional) Lagrange multiplier values are stored
1546  // after those that were created by the bulk elements:
1547  lambda+=node_pt(j)->value(Nbulk_value[j]+k)*psi(j,k);
1548  }
1549  }
1550 
1551  // Get vector of coordinates in bulk element
1552  s_bulk=local_coordinate_in_bulk(s);
1553 
1554  // Get unit normal on bulk shell element
1555  bulk_el_pt->get_normal(s_bulk,shell_normal);
1556 
1557  //Premultiply the weights and the Jacobian
1558  double W = w*J;
1559 
1560  // Here's the actual constraint
1561  double constraint_residual=
1562  shell_normal[0]*Normal_to_clamping_plane[0]+
1563  shell_normal[1]*Normal_to_clamping_plane[1]+
1564  shell_normal[2]*Normal_to_clamping_plane[2];
1565 
1566 
1567  // Assemble residual for Lagrange multiplier:
1568  //-------------------------------------------
1569 
1570  //Loop over the number of nodes
1571  for(unsigned j=0;j<n_node;j++)
1572  {
1573  //Loop over the type of degree of freedom
1574  for(unsigned k=0;k<n_type;k++)
1575  {
1576  // Local eqn number. Recall that the
1577  // (additional) Lagrange multiplier values are stored
1578  // after those that were created by the bulk elements:
1579  local_eqn=nodal_local_eqn(j,Nbulk_value[j]+k);
1580  if (local_eqn>=0)
1581  {
1582  residuals[local_eqn]+=constraint_residual*psi(j,k)*W;
1583  }
1584  }
1585  }
1586 
1587 
1588  // Add Lagrange multiplier contribution to bulk equations
1589  //-------------------------------------------------------
1590 
1591  // Derivatives of constraint equations w.r.t. positional values
1592  // is done by finite differencing
1593  double fd_step=1.0e-8;
1594 
1595  //Loop over the number of nodes
1596  for(unsigned j=0;j<n_node;j++)
1597  {
1598  Node* nod_pt=node_pt(j);
1599 
1600  // Loop over directions
1601  for (unsigned i=0;i<3;i++)
1602  {
1603  //Loop over the type of degree of freedom
1604  for(unsigned k=0;k<4;k++)
1605  {
1606  // Local eqn number: Node, type, direction
1607  local_eqn=position_local_eqn(j,k,i);
1608  if (local_eqn>=0)
1609  {
1610 
1611  // Backup
1612  double backup=nod_pt->x_gen(k,i);
1613 
1614  // Increment
1615  nod_pt->x_gen(k,i)+=fd_step;
1616 
1617  // Re-evaluate unit normal on bulk shell element
1618  bulk_el_pt->get_normal(s_bulk,shell_normal_plus);
1619 
1620  // Re-evaluate the constraint
1621  double constraint_residual_plus=
1622  shell_normal_plus[0]*Normal_to_clamping_plane[0]+
1623  shell_normal_plus[1]*Normal_to_clamping_plane[1]+
1624  shell_normal_plus[2]*Normal_to_clamping_plane[2];
1625 
1626  // Derivarive of constraint w.r.t. to current discrete
1627  // displacement
1628  double dres_dx=(constraint_residual_plus-constraint_residual)/
1629  fd_step;
1630 
1631  // Add to residual
1632  residuals[local_eqn]+=lambda*dres_dx*W;
1633 
1634  // Reset
1635  nod_pt->x_gen(k,i)=backup;
1636 
1637  }
1638  }
1639  }
1640  }
1641  } //End of loop over the integration points
1642 
1643  // Loop over integration points
1644 
1645 }
1646 
1647 
1648 
1649 //===========================================================================
1650 /// Output function
1651 //===========================================================================
1653  output(std::ostream &outfile, const unsigned &n_plot)
1654 {
1655 
1656  // Get bulk element
1657  HermiteShellElement* bulk_el_pt=
1658  dynamic_cast<HermiteShellElement*>(bulk_element_pt());
1659 
1660  //Local coord
1661  Vector<double> s(1);
1662 
1663  //Coordinate in bulk
1664  Vector<double> s_bulk(2);
1665 
1666  // Normal vector
1667  Vector<double> shell_normal(3);
1668 
1669  // # of nodes, # of dof types
1670  Shape psi(2,2);
1671 
1672  //Tecplot header info
1673  outfile << "ZONE I=" << n_plot << std::endl;
1674 
1675  //Loop over plot points
1676  for(unsigned l=0;l<n_plot;l++)
1677  {
1678  s[0] = -1.0 + l*2.0/(n_plot-1);
1679 
1680  // Get vector of coordinates in bulk element
1681  s_bulk=local_coordinate_in_bulk(s);
1682 
1683  // Get unit normal on bulk shell element
1684  bulk_el_pt->get_normal(s_bulk,shell_normal);
1685 
1686  // Get shape function
1687  shape(s,psi);
1688 
1689  // Assemble the Lagrange multiplier
1690  double lambda=0.0;
1691  for(unsigned j=0;j<2;j++)
1692  {
1693  for(unsigned k=0;k<2;k++)
1694  {
1695  // The (additional) Lagrange multiplier values are stored
1696  // after those that were created by the bulk elements:
1697  lambda+=node_pt(j)->value(Nbulk_value[j]+k)*psi(j,k);
1698  }
1699  }
1700 
1701  double dot_product=
1702  shell_normal[0]*Normal_to_clamping_plane[0]+
1703  shell_normal[1]*Normal_to_clamping_plane[1]+
1704  shell_normal[2]*Normal_to_clamping_plane[2];
1705 
1706  //Output stuff
1707  outfile << s[0] << " "
1708  << interpolated_xi(s,0) << " "
1709  << interpolated_xi(s,1) << " "
1710  << interpolated_x(s,0) << " "
1711  << interpolated_x(s,1) << " "
1712  << interpolated_x(s,2) << " "
1713  << J_eulerian(s) << " "
1714  << bulk_el_pt->interpolated_xi(s_bulk,0) << " "
1715  << bulk_el_pt->interpolated_xi(s_bulk,1) << " "
1716  << bulk_el_pt->interpolated_x(s_bulk,0) << " "
1717  << bulk_el_pt->interpolated_x(s_bulk,1) << " "
1718  << bulk_el_pt->interpolated_x(s_bulk,2) << " "
1719  << lambda << " "
1720  << shell_normal[0] << " "
1721  << shell_normal[1] << " "
1722  << shell_normal[2] << " "
1723  << Normal_to_clamping_plane[0] << " "
1724  << Normal_to_clamping_plane[1] << " "
1725  << Normal_to_clamping_plane[2] << " "
1726  << dot_product << " "
1727  << std::endl;
1728  }
1729 
1730 
1731  // Initialise error
1732  const bool check_error=false;
1733  if (check_error)
1734  {
1735  double max_error=0.0;
1736 
1737  //Set # of integration points
1738  const unsigned n_intpt = integral_pt()->nweight();
1739 
1740  //Loop over the integration points
1741  for(unsigned ipt=0;ipt<n_intpt;ipt++)
1742  {
1743 
1744  // Get integration point in local coordinates
1745  s[0]=integral_pt()->knot(ipt,0);
1746 
1747  // Get eulerian jacobian via s
1748  double J_via_s = J_eulerian(s);
1749 
1750  // Get eulerian jacobian via ipt
1751  double J_via_ipt = J_eulerian_at_knot(ipt);
1752 
1753  double error=std::fabs(J_via_s-J_via_ipt);
1754  if (error>max_error) max_error=error;
1755  }
1756  if (max_error>1.0e-14)
1757  {
1758  oomph_info << "Warning: Max. error between J_via_s J_via_ipt "
1759  << max_error << std::endl;
1760  }
1761  }
1762 }
1763 
1764 //=============================================================================
1765 /// Create a list of pairs for all unknowns in this element,
1766 /// so that the first entry in each pair contains the global equation
1767 /// number of the unknown, while the second one contains the number
1768 /// of the "DOF types" that this unknown is associated with.
1769 /// (Function can obviously only be called if the equation numbering
1770 /// scheme has been set up.)
1771 /// This element is only in charge of the solid dofs.
1772 //=============================================================================
1774  std::list<std::pair<unsigned long,unsigned> >& dof_lookup_list) const
1775 {
1776 
1777  // temporary pair (used to store dof lookup prior to being added to list)
1778  std::pair<unsigned long,unsigned> dof_lookup;
1779 
1780  // number of nodes
1781  const unsigned n_node = this->nnode();
1782 
1783  //Get the number of position dofs and dimensions at the node
1784  const unsigned n_position_type = nnodal_position_type();
1785  const unsigned nodal_dim = nodal_dimension();
1786 
1787  //Integer storage for local unknown
1788  int local_unknown=0;
1789 
1790  //Loop over the nodes
1791  for(unsigned n=0;n<n_node;n++)
1792  {
1793 
1794  Node* nod_pt = this->node_pt(n);
1795 
1796  //Loop over position dofs
1797  for(unsigned k=0;k<n_position_type;k++)
1798  {
1799  //Loop over dimension
1800  for(unsigned i=0;i<nodal_dim;i++)
1801  {
1802  //If the variable is free
1803  local_unknown = position_local_eqn(n,k,i);
1804 
1805  // ignore pinned values
1806  if (local_unknown >= 0)
1807  {
1808  // store dof lookup in temporary pair: First entry in pair
1809  // is global equation number; second entry is dof type
1810  dof_lookup.first = this->eqn_number(local_unknown);
1811  dof_lookup.second = 0;
1812 
1813  // add to list
1814  dof_lookup_list.push_front(dof_lookup);
1815 
1816  }
1817  }
1818  }
1819 
1820  //Loop over values
1821  unsigned n_val=nod_pt->nvalue();
1822  for(unsigned j=0;j<n_val;j++)
1823  {
1824  //If the variable is free
1825  local_unknown = nodal_local_eqn(n,j);
1826 
1827  // ignore pinned values
1828  if (local_unknown >= 0)
1829  {
1830  // store dof lookup in temporary pair: First entry in pair
1831  // is global equation number; second entry is dof type
1832  dof_lookup.first = this->eqn_number(local_unknown);
1833  dof_lookup.second = 0;
1834 
1835  // add to list
1836  dof_lookup_list.push_front(dof_lookup);
1837 
1838  }
1839  }
1840  }
1841 }
1842 
1843 
1844 }
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
double load_rate_of_work()
Get integral of instantaneous rate of work done on the wall due to the load returned by the virtual f...
Class of matrices containing doubles, and stored as a DenseMatrix<double>, but with solving functiona...
Definition: matrices.h:1234
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Fill in the element's contribution to its residual vector.
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...
static double Default_lambda_sq_value
Static default value for the timescale ratio (1.0 for natural scaling)
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
void output_with_time_dep_quantities(std::ostream &outfile, const unsigned &n_plot)
Output position veloc and accel.
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 J_eulerian_at_knot(const unsigned &ipt) const
Return the Jacobian of the mapping from local to global coordinates at the ipt-th integration point O...
Definition: elements.cc:5121
unsigned ndim() const
Access function to # of Eulerian coordinates.
Definition: geom_objects.h:185
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 Zero_traction_fct(const Vector< double > &xi, const Vector< double > &x, const Vector< double > &N, Vector< double > &load)
Default load function (zero traction)
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Return the jacobian is calculated by finite differences by default,.
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
const double & h() const
Return the wall thickness to undeformed radius ratio.
virtual double interpolated_dxdt(const Vector< double > &s, const unsigned &i, const unsigned &t)
Return t-th time-derivative of the i-th FE-interpolated Eulerian coordinate at local coordinate s...
Definition: elements.cc:4400
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
double 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:2250
const double & lambda_sq() const
Return the timescale ratio (non-dimensional density)
OomphInfo oomph_info
unsigned nodal_dimension() const
Return the required Eulerian dimension of the nodes in this element.
Definition: elements.h:2358
bool Ignore_membrane_terms
Boolean flag to ignore membrane terms.
e
Definition: cfortran.h:575
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
unsigned nvalue() const
Return number of values stored in data object (incl pinned ones).
Definition: nodes.h:448
void dposition_dlagrangian_at_local_coordinate(const Vector< double > &s, DenseMatrix< double > &drdxi) const
Derivative of position vector w.r.t. the SolidFiniteElement's.
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
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
static double Zero_prestress
Static default for prestress (set to zero)
double calculate_contravariant(double A[2][2], double Aup[2][2])
Invert a DIM by DIM matrix.
GeomObject * Undeformed_midplane_pt
Pointer to the GeomObject that specifies the undeformed midplane of the shell.
std::pair< double, double > get_strain_and_bend(const Vector< double > &s, DenseDoubleMatrix &strain, DenseDoubleMatrix &bend)
Get strain and bending tensors; returns pair comprising the determinant of the undeformed (*...
double d2shape_lagrangian(const Vector< double > &s, Shape &psi, DShape &dpsidxi, DShape &d2psidxi) const
Compute the geometric shape functions and also first and second derivatives w.r.t. Lagrangian coordinates at local coordinate s; Returns Jacobian of mapping from Lagrangian to local coordinates. Numbering: 1D: d2pidxi(i,0) = 2D: d2psidxi(i,0) = d2psidxi(i,1) = d2psidxi(i,2) = 3D: d2psidxi(i,0) = d2psidxi(i,1) = d2psidxi(i,2) = d2psidxi(i,3) = d2psidxi(i,4) = d2psidxi(i,5) = .
Definition: elements.cc:6456
void bulk_position_type_resize(const unsigned &i)
Resize the storage for bulk_position_type to i entries.
Definition: elements.h:4547
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
unsigned nlagrangian() const
Access function to # of Lagrangian coordinates.
Definition: geom_objects.h:182
static double Default_nu_value
Static default value for the Poisson's ratio.
A Rank 3 Tensor class.
Definition: matrices.h:1337
static double Default_h_value
Static default value for the thickness ratio.
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
DenseMatrix< double * > Prestress_pt
Pointer to membrane pre-stress terms – should probably generalise this to function pointers at some p...
static char t char * s
Definition: cfortran.h:572
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. coordinate and normal vector and...
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
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
void shape(const Vector< double > &s, Shape &psi) const
Calculate the geometric shape functions at local coordinate s. Set any "superfluous" shape functions ...
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
double & x_gen(const unsigned &k, const unsigned &i)
Reference to the generalised position x(k,i).
Definition: nodes.h:1055
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
double J_eulerian(const Vector< double > &s) const
Return the Jacobian of mapping from local to global coordinates at local position s...
Definition: elements.cc:5033
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...
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
unsigned & bulk_position_type(const unsigned &i)
Return the position type in the "bulk" element that corresponds to position type i on the FaceElement...
Definition: elements.h:4526
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
double 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:3646
Vector< double > Normal_to_clamping_plane
Normal vector to the clamping plane.
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
double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s. Overloaded to get information from bulk...
Definition: elements.h:4277
unsigned ndof() const
Return the number of equations/dofs in the element.
Definition: elements.h:834
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
double prestress(const unsigned &i, const unsigned &j)
Return (i,j)-th component of second Piola Kirchhoff membrane prestress.
void output(std::ostream &outfile)
Overload the output function.
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
FiniteElement *& bulk_element_pt()
Pointer to higher-dimensional "bulk" element.
Definition: elements.h:4470
void fill_in_contribution_to_residuals_shell(Vector< double > &residuals)
Return the residuals for the equations of KL shell theory.
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
const double & nu() const
Return the Poisson's ratio.
Vector< double > local_coordinate_in_bulk(const Vector< double > &s) const
Return vector of local coordinates in bulk element, given the local coordinates in this FaceElement...
Definition: elements.cc:6032
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
double interpolated_xi(const Vector< double > &s, const unsigned &i) const
Return i-th FE-interpolated Lagrangian coordinate xi[i] at local coordinate s. Overloaded from SolidF...
Definition: elements.h:4638
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
ClampedHermiteShellBoundaryConditionElement()
Broken empty constructor.
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 load_vector_for_rate_of_work_computation(const unsigned &intpt, const Vector< double > &xi, const Vector< double > &x, const Vector< double > &N, Vector< double > &load)
Get the load vector for the computation of the rate of work done by the load. Here we simply forward ...
void get_energy(double &pot_en, double &kin_en)
Get potential (strain) and kinetic energy of the element.
void get_normal(const Vector< double > &s, Vector< double > &N)
Get normal vector.
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