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