hermite_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 the QHermiteElement classes
27
28// oomph-lib header
29#include "hermite_elements.h"
31
32namespace oomph
33{
34 /// /////////////////////////////////////////////////////////////
35 // 1D Hermite elements
36 /// /////////////////////////////////////////////////////////////
37
38 //=======================================================================
39 /// Assign the static Default_integration_scheme
40 //=======================================================================
41 template<unsigned DIM>
43
44 //=======================================================================
45 /// Shape function for specific QHermiteElement<1>
46 //=======================================================================
47 template<>
49 {
50 // Local storage
51 double Psi[2][2];
52 // Call the OneDimensional Shape functions
53 OneDimHermite::shape(s[0], Psi);
54
55 // Loop over the number of nodes
56 for (unsigned l = 0; l < 2; l++)
57 {
58 // Loop over the number of dofs
59 for (unsigned k = 0; k < 2; k++)
60 {
61 psi(l, k) = Psi[l][k];
62 }
63 }
64 }
65
66 //=======================================================================
67 /// Derivatives of shape functions for specific QHermiteElement<1>
68 //=======================================================================
69 template<>
71 Shape& psi,
72 DShape& dpsids) const
73 {
74 // Local storage
75 double Psi[2][2], DPsi[2][2];
76 // Call the OneDimensional Shape functions
77 OneDimHermite::shape(s[0], Psi);
78 OneDimHermite::dshape(s[0], DPsi);
79
80 // Loop over number of nodes
81 for (unsigned l = 0; l < 2; l++)
82 {
83 // Loop over number of dofs
84 for (unsigned k = 0; k < 2; k++)
85 {
86 psi(l, k) = Psi[l][k];
87 dpsids(l, k, 0) = DPsi[l][k];
88 }
89 }
90 }
91
92 //=======================================================================
93 /// Derivatives and second derivatives of shape functions for specific
94 /// QHermiteElement<1>
95 /// d2psids(i,0) = \f$ d^2 \psi_j / d s^2 \f$
96 //=======================================================================
97 template<>
99 Shape& psi,
100 DShape& dpsids,
101 DShape& d2psids) const
102 {
103 // Local storage
104 double Psi[2][2], DPsi[2][2], D2Psi[2][2];
105 // Call the OneDimensional Shape functions
106 OneDimHermite::shape(s[0], Psi);
107 OneDimHermite::dshape(s[0], DPsi);
108 OneDimHermite::d2shape(s[0], D2Psi);
109
110 // Loop over number of nodes
111 for (unsigned l = 0; l < 2; l++)
112 {
113 // Loop over number of dofs
114 for (unsigned k = 0; k < 2; k++)
115 {
116 psi(l, k) = Psi[l][k];
117 dpsids(l, k, 0) = DPsi[l][k];
118 d2psids(l, k, 0) = D2Psi[l][k];
119 }
120 }
121 }
122
123
124 //=======================================================================
125 /// The output function for general 1D QHermiteElements
126 //=======================================================================
127 template<>
128 void QHermiteElement<1>::output(std::ostream& outfile)
129 {
130 // Tecplot header info
131 outfile << "ZONE I=" << 2 << std::endl;
132
133 // Find the dimension of the node
134 unsigned n_dim = this->nodal_dimension();
135
136 // Loop over element nodes
137 for (unsigned l = 0; l < 2; l++)
138 {
139 // Loop over the dimensions and output the position
140 for (unsigned i = 0; i < n_dim; i++)
141 {
142 outfile << node_pt(l)->x(i) << " ";
143 }
144
145 // Find the number of types of dof stored at each node
146 unsigned n_position_type = node_pt(l)->nposition_type();
147 // Loop over the additional positional dofs
148 for (unsigned k = 1; k < n_position_type; k++)
149 {
150 for (unsigned i = 0; i < n_dim; i++)
151 {
152 outfile << node_pt(l)->x_gen(k, i) << " ";
153 }
154 }
155
156 // Find out how many data values at the node
157 unsigned initial_nvalue = node_pt(l)->nvalue();
158 // Lopp over the data and output whether pinned or not
159 for (unsigned i = 0; i < initial_nvalue; i++)
160 {
161 outfile << node_pt(l)->is_pinned(i) << " ";
162 }
163 outfile << std::endl;
164 }
165 outfile << std::endl;
166 }
167
168 //=======================================================================
169 /// The output function for n_plot points in each coordinate direction
170 //=======================================================================
171 template<>
172 void QHermiteElement<1>::output(std::ostream& outfile, const unsigned& n_plot)
173 {
174 // Local variables
175 Vector<double> s(1);
176
177 // Tecplot header info
178 outfile << "ZONE I=" << n_plot << std::endl;
179
180 // Find the dimension of the first node
181 unsigned n_dim = this->nodal_dimension();
182
183 // Loop over plot points
184 for (unsigned l = 0; l < n_plot; l++)
185 {
186 s[0] = -1.0 + l * 2.0 / (n_plot - 1);
187
188 // Output the coordinates
189 for (unsigned i = 0; i < n_dim; i++)
190 {
191 outfile << interpolated_x(s, i) << " ";
192 }
193 outfile << std::endl;
194 }
195 outfile << std::endl;
196 }
197
198
199 //=======================================================================
200 /// The C-style output function for general 1D QHermiteElements
201 //=======================================================================
202 template<>
203 void QHermiteElement<1>::output(FILE* file_pt)
204 {
205 // Tecplot header info
206 fprintf(file_pt, "ZONE I=2\n");
207
208 // Find the dimension of the nodes
209 unsigned n_dim = this->nodal_dimension();
210
211 // Loop over element nodes
212 for (unsigned l = 0; l < 2; l++)
213 {
214 // Loop over the dimensions and output the position
215 for (unsigned i = 0; i < n_dim; i++)
216 {
217 fprintf(file_pt, "%g ", node_pt(l)->x(i));
218 }
219
220 // Find the number of types of dof stored at each node
221 unsigned n_position_type = node_pt(l)->nposition_type();
222 // Loop over the additional positional dofs
223 for (unsigned k = 1; k < n_position_type; k++)
224 {
225 for (unsigned i = 0; i < n_dim; i++)
226 {
227 fprintf(file_pt, "%g ", node_pt(l)->x_gen(k, i));
228 }
229 }
230
231 // Find out how many data values at the node
232 unsigned initial_nvalue = node_pt(l)->nvalue();
233 // Lopp over the data and output whether pinned or not
234 for (unsigned i = 0; i < initial_nvalue; i++)
235 {
236 fprintf(file_pt, "%i ", node_pt(l)->is_pinned(i));
237 }
238 fprintf(file_pt, "\n");
239 }
240 fprintf(file_pt, "\n");
241 }
242
243
244 //=======================================================================
245 /// The C-style output function for n_plot points in each coordinate direction
246 //=======================================================================
247 template<>
248 void QHermiteElement<1>::output(FILE* file_pt, const unsigned& n_plot)
249 {
250 // Local variables
251 Vector<double> s(1);
252
253 // Tecplot header info
254 fprintf(file_pt, "ZONE I=%i \n", n_plot);
255
256 // Find the dimension of the first node
257 unsigned n_dim = this->nodal_dimension();
258
259 // Loop over element nodes
260 for (unsigned l = 0; l < n_plot; l++)
261 {
262 s[0] = -1.0 + l * 2.0 / (n_plot - 1);
263 // Output the coordinates
264 for (unsigned i = 0; i < n_dim; i++)
265 {
266 fprintf(file_pt, "%g ", interpolated_x(s, i));
267 }
268 fprintf(file_pt, "\n");
269 }
270 }
271
272
273 //===========================================================
274 /// Function to setup geometrical information for lower-dimensional
275 /// FaceElements (single node elements)
276 //===========================================================
277 template<>
278 void QHermiteElement<1>::build_face_element(const int& face_index,
279 FaceElement* face_element_pt)
280 {
281 // Set the nodal dimension from the "bulk"
282 face_element_pt->set_nodal_dimension(node_pt(0)->ndim());
283
284 // Set the pointer to the "bulk" element
285 face_element_pt->bulk_element_pt() = this;
286
287#ifdef OOMPH_HAS_MPI
288 // Pass on non-halo ID
289 face_element_pt->set_halo(Non_halo_proc_ID);
290#endif
291
292 // Resize the storage for the original number of values at the (one and
293 // only) node of the face element.
294 face_element_pt->nbulk_value_resize(1);
295
296 // Resize the storage for the bulk node number corresponding to the (one
297 // and only) node of the face element
298 face_element_pt->bulk_node_number_resize(1);
299
300 // Set the face index in the face element
301 face_element_pt->face_index() = face_index;
302
303 // Now set up the node pointer
304 // The convention is that the "normal", the coordinate, should always point
305 // out of the element
306 switch (face_index)
307 {
308 // Bottom, normal sign is negative (coordinate points into element)
309 case (-1):
310 face_element_pt->node_pt(0) = node_pt(0);
311 face_element_pt->bulk_node_number(0) = 0;
312 face_element_pt->normal_sign() = -1;
313
314 // Set the pointer to the function that determines the bulk coordinates
315 // in the face element
316 face_element_pt->face_to_bulk_coordinate_fct_pt() =
318
319 // Set the pointer to the function that determines the mapping of
320 // derivatives
321 face_element_pt->bulk_coordinate_derivatives_fct_pt() =
323
324 // Set the number of values stored when the node is part of the "bulk"
325 // element.
326 face_element_pt->nbulk_value(0) = required_nvalue(0);
327 break;
328
329 // Top, normal sign is positive (coordinate points out of element)
330 case (1):
331 face_element_pt->node_pt(0) = node_pt(1);
332 face_element_pt->bulk_node_number(0) = 1;
333 face_element_pt->normal_sign() = +1;
334
335 // Set the pointer to the function that determines the bulk coordinates
336 // in the face element
337 face_element_pt->face_to_bulk_coordinate_fct_pt() =
339
340 // Set the pointer to the function that determines the mapping of
341 // derivatives
342 face_element_pt->bulk_coordinate_derivatives_fct_pt() =
344
345 // Set the number of values stored when the node is part of the "bulk"
346 // element.
347 face_element_pt->nbulk_value(0) = required_nvalue(1);
348 break;
349
350 // Other cases
351 default:
352 std::ostringstream error_message;
353 error_message << "Face_index should only take "
354 << "the values +/-1, not " << face_index << std::endl;
355
356 throw OomphLibError(error_message.str(),
357 OOMPH_CURRENT_FUNCTION,
358 OOMPH_EXCEPTION_LOCATION);
359 }
360 }
361
362 /// /////////////////////////////////////////////////////////////
363 // 2D Hermite elements
364 /// /////////////////////////////////////////////////////////////
365
366
367 //=======================================================================
368 /// Shape function for specific QHermiteElement<2>
369 //=======================================================================
370 template<>
372 {
373 // Local storage
374 double Psi[2][2][2];
375
376 // Call the OneDimensional Shape functions
377 OneDimHermite::shape(s[0], Psi[0]);
378 OneDimHermite::shape(s[1], Psi[1]);
379
380 // Set up the two dimensional shape functions
381 // Set up the functions at corner 0
382 // psi_0 = 1 when s1 = 0, s2 = 0
383 psi(0, 0) = Psi[0][0][0] * Psi[1][0][0];
384 // dpsi_0/ds1 = 1 when s1 = 0, s2 = 0
385 psi(0, 1) = Psi[0][0][1] * Psi[1][0][0];
386 // dpsi_0/ds2 = 1 when s1 = 0, s2 = 0
387 psi(0, 2) = Psi[0][0][0] * Psi[1][0][1];
388 // dpsi_0/ds2ds1 = 1 when s1 = 0, s2 = 0
389 psi(0, 3) = Psi[0][0][1] * Psi[1][0][1];
390
391 // Set up the functions at corner 1
392 // psi_1 = 1 when s1 = 1, s2 = 0
393 psi(1, 0) = Psi[0][1][0] * Psi[1][0][0];
394 // dpsi_1/ds1 = 1 when s1 = 1, s2 = 0
395 psi(1, 1) = Psi[0][1][1] * Psi[1][0][0];
396 // dpsi_1/ds2 = 1 when s1 = 1, s2 = 0
397 psi(1, 2) = Psi[0][1][0] * Psi[1][0][1];
398 // dpsi_1/ds1ds2 = 1 when s1 = 1, s2 = 0
399 psi(1, 3) = Psi[0][1][1] * Psi[1][0][1];
400
401 // Set up the functions at the corner 2
402 // psi_2 = 1 when s1 = 0, s2 = 1
403 psi(2, 0) = Psi[0][0][0] * Psi[1][1][0];
404 // dpsi_2/ds1 = 1 when s1 = 0, s2 = 1
405 psi(2, 1) = Psi[0][0][1] * Psi[1][1][0];
406 // dpsi_2/ds2 = 1 when s1 = 0, s2 = 1
407 psi(2, 2) = Psi[0][0][0] * Psi[1][1][1];
408 // dpsi_2/ds2ds1 = 1 when s1 = 0, s2 = 1
409 psi(2, 3) = Psi[0][0][1] * Psi[1][1][1];
410
411 // Set up the functions at corner 3
412 // psi_3 = 1 when s1 = 1, s2 = 1
413 psi(3, 0) = Psi[0][1][0] * Psi[1][1][0];
414 // dpsi_3/ds1 = 1 when s1 = 1, s2 = 1
415 psi(3, 1) = Psi[0][1][1] * Psi[1][1][0];
416 // dpsi_3/ds2 = 1 when s1 = 1, s2 = 1
417 psi(3, 2) = Psi[0][1][0] * Psi[1][1][1];
418 // dpsi_3/ds1ds2 = 1 when s1 = 1, s2 = 1
419 psi(3, 3) = Psi[0][1][1] * Psi[1][1][1];
420 }
421
422 //=======================================================================
423 /// Derivatives of shape functions for specific QHermiteElement<2>
424 //=======================================================================
425 template<>
427 Shape& psi,
428 DShape& dpsids) const
429 {
430 // Local storage
431 double Psi[2][2][2];
432 double DPsi[2][2][2];
433 // Call the OneDimensional Shape functions
434 OneDimHermite::shape(s[0], Psi[0]);
435 OneDimHermite::shape(s[1], Psi[1]);
436 OneDimHermite::dshape(s[0], DPsi[0]);
437 OneDimHermite::dshape(s[1], DPsi[1]);
438
439
440 // Set up the two dimensional shape functions
441 // Set up the functions at corner 0
442 // psi_0 = 1 when s1 = 0, s2 = 0
443 psi(0, 0) = Psi[0][0][0] * Psi[1][0][0];
444 // dpsi_0/ds1 = 1 when s1 = 0, s2 = 0
445 psi(0, 1) = Psi[0][0][1] * Psi[1][0][0];
446 // dpsi_0/ds2 = 1 when s1 = 0, s2 = 0
447 psi(0, 2) = Psi[0][0][0] * Psi[1][0][1];
448 // dpsi_0/ds2ds1 = 1 when s1 = 0, s2 = 0
449 psi(0, 3) = Psi[0][0][1] * Psi[1][0][1];
450
451 // Set up the functions at corner 1
452 // psi_1 = 1 when s1 = 1, s2 = 0
453 psi(1, 0) = Psi[0][1][0] * Psi[1][0][0];
454 // dpsi_1/ds1 = 1 when s1 = 1, s2 = 0
455 psi(1, 1) = Psi[0][1][1] * Psi[1][0][0];
456 // dpsi_1/ds2 = 1 when s1 = 1, s2 = 0
457 psi(1, 2) = Psi[0][1][0] * Psi[1][0][1];
458 // dpsi_1/ds1ds2 = 1 when s1 = 1, s2 = 0
459 psi(1, 3) = Psi[0][1][1] * Psi[1][0][1];
460
461 // Set up the functions at the corner 2
462 // psi_2 = 1 when s1 = 0, s2 = 1
463 psi(2, 0) = Psi[0][0][0] * Psi[1][1][0];
464 // dpsi_2/ds1 = 1 when s1 = 0, s2 = 1
465 psi(2, 1) = Psi[0][0][1] * Psi[1][1][0];
466 // dpsi_2/ds2 = 1 when s1 = 0, s2 = 1
467 psi(2, 2) = Psi[0][0][0] * Psi[1][1][1];
468 // dpsi_2/ds2ds1 = 1 when s1 = 0, s2 = 1
469 psi(2, 3) = Psi[0][0][1] * Psi[1][1][1];
470
471 // Set up the functions at corner 3
472 // psi_3 = 1 when s1 = 1, s2 = 1
473 psi(3, 0) = Psi[0][1][0] * Psi[1][1][0];
474 // dpsi_3/ds1 = 1 when s1 = 1, s2 = 1
475 psi(3, 1) = Psi[0][1][1] * Psi[1][1][0];
476 // dpsi_3/ds2 = 1 when s1 = 1, s2 = 1
477 psi(3, 2) = Psi[0][1][0] * Psi[1][1][1];
478 // dpsi_3/ds1ds2 = 1 when s1 = 1, s2 = 1
479 psi(3, 3) = Psi[0][1][1] * Psi[1][1][1];
480
481 // FIRST DERIVATIVES
482
483 // D/Ds[0]
484
485 // Set up the functions at corner 0
486 dpsids(0, 0, 0) = DPsi[0][0][0] * Psi[1][0][0];
487 dpsids(0, 1, 0) = DPsi[0][0][1] * Psi[1][0][0];
488 dpsids(0, 2, 0) = DPsi[0][0][0] * Psi[1][0][1];
489 dpsids(0, 3, 0) = DPsi[0][0][1] * Psi[1][0][1];
490
491 // Set up the functions at corner 1
492 dpsids(1, 0, 0) = DPsi[0][1][0] * Psi[1][0][0];
493 dpsids(1, 1, 0) = DPsi[0][1][1] * Psi[1][0][0];
494 dpsids(1, 2, 0) = DPsi[0][1][0] * Psi[1][0][1];
495 dpsids(1, 3, 0) = DPsi[0][1][1] * Psi[1][0][1];
496
497 // Set up the functions at the corner 2
498 dpsids(2, 0, 0) = DPsi[0][0][0] * Psi[1][1][0];
499 dpsids(2, 1, 0) = DPsi[0][0][1] * Psi[1][1][0];
500 dpsids(2, 2, 0) = DPsi[0][0][0] * Psi[1][1][1];
501 dpsids(2, 3, 0) = DPsi[0][0][1] * Psi[1][1][1];
502
503 // Set up the functions at corner 3
504 dpsids(3, 0, 0) = DPsi[0][1][0] * Psi[1][1][0];
505 dpsids(3, 1, 0) = DPsi[0][1][1] * Psi[1][1][0];
506 dpsids(3, 2, 0) = DPsi[0][1][0] * Psi[1][1][1];
507 dpsids(3, 3, 0) = DPsi[0][1][1] * Psi[1][1][1];
508
509 // D/Ds[1]
510
511 // Set up the functions at corner 0
512 dpsids(0, 0, 1) = Psi[0][0][0] * DPsi[1][0][0];
513 dpsids(0, 1, 1) = Psi[0][0][1] * DPsi[1][0][0];
514 dpsids(0, 2, 1) = Psi[0][0][0] * DPsi[1][0][1];
515 dpsids(0, 3, 1) = Psi[0][0][1] * DPsi[1][0][1];
516
517 // Set up the functions at corner 1
518 dpsids(1, 0, 1) = Psi[0][1][0] * DPsi[1][0][0];
519 dpsids(1, 1, 1) = Psi[0][1][1] * DPsi[1][0][0];
520 dpsids(1, 2, 1) = Psi[0][1][0] * DPsi[1][0][1];
521 dpsids(1, 3, 1) = Psi[0][1][1] * DPsi[1][0][1];
522
523 // Set up the functions at the corner 2
524 dpsids(2, 0, 1) = Psi[0][0][0] * DPsi[1][1][0];
525 dpsids(2, 1, 1) = Psi[0][0][1] * DPsi[1][1][0];
526 dpsids(2, 2, 1) = Psi[0][0][0] * DPsi[1][1][1];
527 dpsids(2, 3, 1) = Psi[0][0][1] * DPsi[1][1][1];
528
529 // Set up the functions at corner 3
530 dpsids(3, 0, 1) = Psi[0][1][0] * DPsi[1][1][0];
531 dpsids(3, 1, 1) = Psi[0][1][1] * DPsi[1][1][0];
532 dpsids(3, 2, 1) = Psi[0][1][0] * DPsi[1][1][1];
533 dpsids(3, 3, 1) = Psi[0][1][1] * DPsi[1][1][1];
534 }
535
536 //======================================================================
537 /// Second derivatives of the shape functions wrt local coordinates.
538 /// d2psids(i,0) = \f$ \partial^2 \psi_j / \partial s_0^2 \f$
539 /// d2psids(i,1) = \f$ \partial^2 \psi_j / \partial s_1^2 \f$
540 /// d2psids(i,2) = \f$ \partial^2 \psi_j / \partial s_0 \partial s_1 \f$
541 //======================================================================
542 template<>
544 Shape& psi,
545 DShape& dpsids,
546 DShape& d2psids) const
547 {
548 // Local storage
549 double Psi[2][2][2];
550 double DPsi[2][2][2];
551 double D2Psi[2][2][2];
552
553 // Call the OneDimensional Shape functions
554 OneDimHermite::shape(s[0], Psi[0]);
555 OneDimHermite::shape(s[1], Psi[1]);
557 OneDimHermite::dshape(s[1], DPsi[1]);
558 OneDimHermite::d2shape(s[0], D2Psi[0]);
559 OneDimHermite::d2shape(s[1], D2Psi[1]);
560
561 // Set up the two dimensional shape functions
562 // Set up the functions at corner 0
563 // psi_0 = 1 when s1 = 0, s2 = 0
564 psi(0, 0) = Psi[0][0][0] * Psi[1][0][0];
565 // dpsi_0/ds1 = 1 when s1 = 0, s2 = 0
566 psi(0, 1) = Psi[0][0][1] * Psi[1][0][0];
567 // dpsi_0/ds2 = 1 when s1 = 0, s2 = 0
568 psi(0, 2) = Psi[0][0][0] * Psi[1][0][1];
569 // dpsi_0/ds2ds1 = 1 when s1 = 0, s2 = 0
570 psi(0, 3) = Psi[0][0][1] * Psi[1][0][1];
571
572 // Set up the functions at corner 1
573 // psi_1 = 1 when s1 = 1, s2 = 0
574 psi(1, 0) = Psi[0][1][0] * Psi[1][0][0];
575 // dpsi_1/ds1 = 1 when s1 = 1, s2 = 0
576 psi(1, 1) = Psi[0][1][1] * Psi[1][0][0];
577 // dpsi_1/ds2 = 1 when s1 = 1, s2 = 0
578 psi(1, 2) = Psi[0][1][0] * Psi[1][0][1];
579 // dpsi_1/ds1ds2 = 1 when s1 = 1, s2 = 0
580 psi(1, 3) = Psi[0][1][1] * Psi[1][0][1];
581
582 // Set up the functions at the corner 2
583 // psi_2 = 1 when s1 = 0, s2 = 1
584 psi(2, 0) = Psi[0][0][0] * Psi[1][1][0];
585 // dpsi_2/ds1 = 1 when s1 = 0, s2 = 1
586 psi(2, 1) = Psi[0][0][1] * Psi[1][1][0];
587 // dpsi_2/ds2 = 1 when s1 = 0, s2 = 1
588 psi(2, 2) = Psi[0][0][0] * Psi[1][1][1];
589 // dpsi_2/ds2ds1 = 1 when s1 = 0, s2 = 1
590 psi(2, 3) = Psi[0][0][1] * Psi[1][1][1];
591
592 // Set up the functions at corner 3
593 // psi_3 = 1 when s1 = 1, s2 = 1
594 psi(3, 0) = Psi[0][1][0] * Psi[1][1][0];
595 // dpsi_3/ds1 = 1 when s1 = 1, s2 = 1
596 psi(3, 1) = Psi[0][1][1] * Psi[1][1][0];
597 // dpsi_3/ds2 = 1 when s1 = 1, s2 = 1
598 psi(3, 2) = Psi[0][1][0] * Psi[1][1][1];
599 // dpsi_3/ds1ds2 = 1 when s1 = 1, s2 = 1
600 psi(3, 3) = Psi[0][1][1] * Psi[1][1][1];
601
602 // FIRST DERIVATIVES
603
604 // D/Ds[0]
605
606 // Set up the functions at corner 0
607 dpsids(0, 0, 0) = DPsi[0][0][0] * Psi[1][0][0];
608 dpsids(0, 1, 0) = DPsi[0][0][1] * Psi[1][0][0];
609 dpsids(0, 2, 0) = DPsi[0][0][0] * Psi[1][0][1];
610 dpsids(0, 3, 0) = DPsi[0][0][1] * Psi[1][0][1];
611
612 // Set up the functions at corner 1
613 dpsids(1, 0, 0) = DPsi[0][1][0] * Psi[1][0][0];
614 dpsids(1, 1, 0) = DPsi[0][1][1] * Psi[1][0][0];
615 dpsids(1, 2, 0) = DPsi[0][1][0] * Psi[1][0][1];
616 dpsids(1, 3, 0) = DPsi[0][1][1] * Psi[1][0][1];
617
618 // Set up the functions at the corner 2
619 dpsids(2, 0, 0) = DPsi[0][0][0] * Psi[1][1][0];
620 dpsids(2, 1, 0) = DPsi[0][0][1] * Psi[1][1][0];
621 dpsids(2, 2, 0) = DPsi[0][0][0] * Psi[1][1][1];
622 dpsids(2, 3, 0) = DPsi[0][0][1] * Psi[1][1][1];
623
624 // Set up the functions at corner 3
625 dpsids(3, 0, 0) = DPsi[0][1][0] * Psi[1][1][0];
626 dpsids(3, 1, 0) = DPsi[0][1][1] * Psi[1][1][0];
627 dpsids(3, 2, 0) = DPsi[0][1][0] * Psi[1][1][1];
628 dpsids(3, 3, 0) = DPsi[0][1][1] * Psi[1][1][1];
629
630 // D/Ds[1]
631
632 // Set up the functions at corner 0
633 dpsids(0, 0, 1) = Psi[0][0][0] * DPsi[1][0][0];
634 dpsids(0, 1, 1) = Psi[0][0][1] * DPsi[1][0][0];
635 dpsids(0, 2, 1) = Psi[0][0][0] * DPsi[1][0][1];
636 dpsids(0, 3, 1) = Psi[0][0][1] * DPsi[1][0][1];
637
638 // Set up the functions at corner 1
639 dpsids(1, 0, 1) = Psi[0][1][0] * DPsi[1][0][0];
640 dpsids(1, 1, 1) = Psi[0][1][1] * DPsi[1][0][0];
641 dpsids(1, 2, 1) = Psi[0][1][0] * DPsi[1][0][1];
642 dpsids(1, 3, 1) = Psi[0][1][1] * DPsi[1][0][1];
643
644 // Set up the functions at the corner 2
645 dpsids(2, 0, 1) = Psi[0][0][0] * DPsi[1][1][0];
646 dpsids(2, 1, 1) = Psi[0][0][1] * DPsi[1][1][0];
647 dpsids(2, 2, 1) = Psi[0][0][0] * DPsi[1][1][1];
648 dpsids(2, 3, 1) = Psi[0][0][1] * DPsi[1][1][1];
649
650 // Set up the functions at corner 3
651 dpsids(3, 0, 1) = Psi[0][1][0] * DPsi[1][1][0];
652 dpsids(3, 1, 1) = Psi[0][1][1] * DPsi[1][1][0];
653 dpsids(3, 2, 1) = Psi[0][1][0] * DPsi[1][1][1];
654 dpsids(3, 3, 1) = Psi[0][1][1] * DPsi[1][1][1];
655
656 // SECOND DERIVATIVES
657 // Convention: index 0 is d^2/ds[0]^2,
658 // index 1 is d^2/ds[1]^2,
659 // index 2 is the mixed derivative
660
661 // D^2/Ds[0]^2
662
663 // Set up the functions at corner 0
664 d2psids(0, 0, 0) = D2Psi[0][0][0] * Psi[1][0][0];
665 d2psids(0, 1, 0) = D2Psi[0][0][1] * Psi[1][0][0];
666 d2psids(0, 2, 0) = D2Psi[0][0][0] * Psi[1][0][1];
667 d2psids(0, 3, 0) = D2Psi[0][0][1] * Psi[1][0][1];
668
669 // Set up the functions at corner 1
670 d2psids(1, 0, 0) = D2Psi[0][1][0] * Psi[1][0][0];
671 d2psids(1, 1, 0) = D2Psi[0][1][1] * Psi[1][0][0];
672 d2psids(1, 2, 0) = D2Psi[0][1][0] * Psi[1][0][1];
673 d2psids(1, 3, 0) = D2Psi[0][1][1] * Psi[1][0][1];
674
675 // Set up the functions at the corner 2
676 d2psids(2, 0, 0) = D2Psi[0][0][0] * Psi[1][1][0];
677 d2psids(2, 1, 0) = D2Psi[0][0][1] * Psi[1][1][0];
678 d2psids(2, 2, 0) = D2Psi[0][0][0] * Psi[1][1][1];
679 d2psids(2, 3, 0) = D2Psi[0][0][1] * Psi[1][1][1];
680
681 // Set up the functions at corner 3
682 d2psids(3, 0, 0) = D2Psi[0][1][0] * Psi[1][1][0];
683 d2psids(3, 1, 0) = D2Psi[0][1][1] * Psi[1][1][0];
684 d2psids(3, 2, 0) = D2Psi[0][1][0] * Psi[1][1][1];
685 d2psids(3, 3, 0) = D2Psi[0][1][1] * Psi[1][1][1];
686
687 // D^2/Ds[1]^2
688
689 // Set up the functions at corner 0
690 d2psids(0, 0, 1) = Psi[0][0][0] * D2Psi[1][0][0];
691 d2psids(0, 1, 1) = Psi[0][0][1] * D2Psi[1][0][0];
692 d2psids(0, 2, 1) = Psi[0][0][0] * D2Psi[1][0][1];
693 d2psids(0, 3, 1) = Psi[0][0][1] * D2Psi[1][0][1];
694
695 // Set up the functions at corner 1
696 d2psids(1, 0, 1) = Psi[0][1][0] * D2Psi[1][0][0];
697 d2psids(1, 1, 1) = Psi[0][1][1] * D2Psi[1][0][0];
698 d2psids(1, 2, 1) = Psi[0][1][0] * D2Psi[1][0][1];
699 d2psids(1, 3, 1) = Psi[0][1][1] * D2Psi[1][0][1];
700
701 // Set up the functions at the corner 2
702 d2psids(2, 0, 1) = Psi[0][0][0] * D2Psi[1][1][0];
703 d2psids(2, 1, 1) = Psi[0][0][1] * D2Psi[1][1][0];
704 d2psids(2, 2, 1) = Psi[0][0][0] * D2Psi[1][1][1];
705 d2psids(2, 3, 1) = Psi[0][0][1] * D2Psi[1][1][1];
706
707 // Set up the functions at corner 3
708 d2psids(3, 0, 1) = Psi[0][1][0] * D2Psi[1][1][0];
709 d2psids(3, 1, 1) = Psi[0][1][1] * D2Psi[1][1][0];
710 d2psids(3, 2, 1) = Psi[0][1][0] * D2Psi[1][1][1];
711 d2psids(3, 3, 1) = Psi[0][1][1] * D2Psi[1][1][1];
712
713 // D^2/Ds[0]Ds[1]
714
715 // Set up the functions at corner 0
716 d2psids(0, 0, 2) = DPsi[0][0][0] * DPsi[1][0][0];
717 d2psids(0, 1, 2) = DPsi[0][0][1] * DPsi[1][0][0];
718 d2psids(0, 2, 2) = DPsi[0][0][0] * DPsi[1][0][1];
719 d2psids(0, 3, 2) = DPsi[0][0][1] * DPsi[1][0][1];
720
721 // Set up the functions at corner 1
722 d2psids(1, 0, 2) = DPsi[0][1][0] * DPsi[1][0][0];
723 d2psids(1, 1, 2) = DPsi[0][1][1] * DPsi[1][0][0];
724 d2psids(1, 2, 2) = DPsi[0][1][0] * DPsi[1][0][1];
725 d2psids(1, 3, 2) = DPsi[0][1][1] * DPsi[1][0][1];
726
727 // Set up the functions at the corner 2
728 d2psids(2, 0, 2) = DPsi[0][0][0] * DPsi[1][1][0];
729 d2psids(2, 1, 2) = DPsi[0][0][1] * DPsi[1][1][0];
730 d2psids(2, 2, 2) = DPsi[0][0][0] * DPsi[1][1][1];
731 d2psids(2, 3, 2) = DPsi[0][0][1] * DPsi[1][1][1];
732
733 // Set up the functions at corner 3
734 d2psids(3, 0, 2) = DPsi[0][1][0] * DPsi[1][1][0];
735 d2psids(3, 1, 2) = DPsi[0][1][1] * DPsi[1][1][0];
736 d2psids(3, 2, 2) = DPsi[0][1][0] * DPsi[1][1][1];
737 d2psids(3, 3, 2) = DPsi[0][1][1] * DPsi[1][1][1];
738 }
739
740
741 //===========================================================
742 /// The output function for QHermiteElement<2,ORDER>
743 //===========================================================
744 template<>
745 void QHermiteElement<2>::output(std::ostream& outfile)
746 {
747 // Tecplot header info
748 outfile << "ZONE I=" << 2 << ", J=" << 2 << std::endl;
749
750 // Find the dimension of the node
751 unsigned n_dim = this->nodal_dimension();
752
753 // Loop over element nodes
754 for (unsigned l2 = 0; l2 < 2; l2++)
755 {
756 for (unsigned l1 = 0; l1 < 2; l1++)
757 {
758 unsigned l = l2 * 2 + l1;
759
760 // Loop over the dimensions and output the position
761 for (unsigned i = 0; i < n_dim; i++)
762 {
763 outfile << node_pt(l)->x(i) << " ";
764 }
765
766 // Find out number of types of dof stored at each node
767 unsigned n_position_type = node_pt(l)->nposition_type();
768 // Loop over the additional positional dofs
769 for (unsigned k = 1; k < n_position_type; k++)
770 {
771 for (unsigned i = 0; i < n_dim; i++)
772 {
773 outfile << node_pt(l)->x_gen(k, i) << " ";
774 }
775 }
776
777 // Find out how many data values at the node
778 unsigned initial_nvalue = node_pt(l)->nvalue();
779 // Loop over the data and output whether pinned or not
780 for (unsigned i = 0; i < initial_nvalue; i++)
781 {
782 outfile << node_pt(l)->is_pinned(i) << " ";
783 }
784 outfile << std::endl;
785 }
786 }
787 outfile << std::endl;
788 }
789
790 //=======================================================================
791 /// The output function for n_plot points in each coordinate direction
792 //=======================================================================
793 template<>
794 void QHermiteElement<2>::output(std::ostream& outfile, const unsigned& n_plot)
795 {
796 // Local variables
797 Vector<double> s(2);
798
799 // Tecplot header info
800 outfile << "ZONE I=" << n_plot << ", J=" << n_plot << std::endl;
801
802 // Find the dimension of the first node
803 unsigned n_dim = this->nodal_dimension();
804
805 // Loop over plot points
806 for (unsigned l2 = 0; l2 < n_plot; l2++)
807 {
808 s[1] = -1.0 + l2 * 2.0 / (n_plot - 1);
809
810 for (unsigned l1 = 0; l1 < n_plot; l1++)
811 {
812 s[0] = -1.0 + l1 * 2.0 / (n_plot - 1);
813
814 // Output the coordinates
815 for (unsigned i = 0; i < n_dim; i++)
816 {
817 outfile << interpolated_x(s, i) << " ";
818 }
819 }
820 }
821 outfile << std::endl;
822 }
823
824
825 //===========================================================
826 /// The C-style output function for QHermiteElement<2,ORDER>
827 //===========================================================
828 template<>
829 void QHermiteElement<2>::output(FILE* file_pt)
830 {
831 // Tecplot header info
832 fprintf(file_pt, "ZONE I=2, J=2");
833
834 // Find the dimension of the node
835 unsigned n_dim = this->nodal_dimension();
836
837 // Loop over element nodes
838 for (unsigned l2 = 0; l2 < 2; l2++)
839 {
840 for (unsigned l1 = 0; l1 < 2; l1++)
841 {
842 unsigned l = l2 * 2 + l1;
843
844 // Loop over the dimensions and output the position
845 for (unsigned i = 0; i < n_dim; i++)
846 {
847 fprintf(file_pt, "%g ", node_pt(l)->x(i));
848 }
849
850 // Find out number of types of dof stored at each node
851 unsigned n_position_type = node_pt(l)->nposition_type();
852 // Loop over the additional positional dofs
853 for (unsigned k = 1; k < n_position_type; k++)
854 {
855 for (unsigned i = 0; i < n_dim; i++)
856 {
857 fprintf(file_pt, "%g ", node_pt(l)->x_gen(k, i));
858 }
859 }
860
861 // Find out how many data values at the node
862 unsigned initial_nvalue = node_pt(l)->nvalue();
863 // Loop over the data and output whether pinned or not
864 for (unsigned i = 0; i < initial_nvalue; i++)
865 {
866 fprintf(file_pt, "%i ", node_pt(l)->is_pinned(i));
867 }
868 fprintf(file_pt, "\n");
869 }
870 }
871 fprintf(file_pt, "\n");
872 }
873
874 //=======================================================================
875 /// The C-style output function for n_plot points in each coordinate direction
876 //=======================================================================
877 template<>
878 void QHermiteElement<2>::output(FILE* file_pt, const unsigned& n_plot)
879 {
880 // Local variables
881 Vector<double> s(2);
882
883 // Tecplot header info
884 fprintf(file_pt, "ZONE I=%i, J=%i \n", n_plot, n_plot);
885
886 // Find the dimension of the nodes
887 unsigned n_dim = this->nodal_dimension();
888
889 // Loop over plot points
890 for (unsigned l2 = 0; l2 < n_plot; l2++)
891 {
892 s[1] = -1.0 + l2 * 2.0 / (n_plot - 1);
893 for (unsigned l1 = 0; l1 < n_plot; l1++)
894 {
895 s[0] = -1.0 + l1 * 2.0 / (n_plot - 1);
896
897 // Output the coordinates
898 for (unsigned i = 0; i < n_dim; i++)
899 {
900 fprintf(file_pt, "%g ", interpolated_x(s, i));
901 }
902 }
903 }
904 fprintf(file_pt, "\n");
905 }
906
907
908 //=======================================================================
909 /// Function to setup geometrical information for lower-dimensional
910 /// FaceElements (of type QHermiteElement<1>).
911 //=======================================================================
912 template<>
913 void QHermiteElement<2>::build_face_element(const int& face_index,
914 FaceElement* face_element_pt)
915 {
916 // Set the nodal dimension from the "bulk"
917 face_element_pt->set_nodal_dimension(node_pt(0)->ndim());
918
919 // Set the pointer to the "bulk" element
920 face_element_pt->bulk_element_pt() = this;
921
922#ifdef OOMPH_HAS_MPI
923 // Pass on non-halo proc ID
924 face_element_pt->set_halo(Non_halo_proc_ID);
925#endif
926
927 // Resize the bulk_position_type translation scheme to the number of
928 // position types in the 1D element: 2 (position and slope)
929 face_element_pt->bulk_position_type_resize(2);
930
931 // Resize the storage for the original number of values at
932 // the two nodes of the FaceElement
933 face_element_pt->nbulk_value_resize(2);
934
935 // Resize the storage for the bulk node numbers coressponding to the
936 // two nodes of the FaceElement
937 face_element_pt->bulk_node_number_resize(2);
938
939 // Set the face index in the face element
940 // The faces are
941 // +1 East
942 // -1 West
943 // +2 North
944 // -2 South
945
946 // Set the face index in the face element
947 face_element_pt->face_index() = face_index;
948
949 // Now set up the node pointers
950 // The convention here is that interior_tangent X tangent X tangent
951 // is the OUTWARD normal
952 // IMPORTANT NOTE: Need to ensure that numbering is consistent here
953 // i.e. node numbers increase in positive x and y directions as they should
954 // If not, normals will be inward.
955 switch (face_index)
956 {
957 unsigned bulk_number;
958 // West face, normal sign is positive
959 case (-1):
960 // Set the pointer to the bulk coordinate translation scheme
961 face_element_pt->face_to_bulk_coordinate_fct_pt() =
963
964 // Set the pointer to the derivative mappings
965 face_element_pt->bulk_coordinate_derivatives_fct_pt() =
967
968 for (unsigned i = 0; i < 2; i++)
969 {
970 bulk_number = i * 2;
971 face_element_pt->node_pt(i) = node_pt(bulk_number);
972 face_element_pt->bulk_node_number(i) = bulk_number;
973 face_element_pt->normal_sign() = 1;
974 // Set the number of values originally stored at this node
975 face_element_pt->nbulk_value(i) = required_nvalue(bulk_number);
976 // Set the position type for the slope, which is in the s[1]
977 // direction, so is position_type 2 in the "bulk" element
978 face_element_pt->bulk_position_type(1) = 2;
979 }
980 break;
981 // South face, normal sign is positive
982 case (-2):
983 // Set the pointer to the bulk coordinate translation scheme
984 face_element_pt->face_to_bulk_coordinate_fct_pt() =
986
987 // Set the pointer to the derivative mappings
988 face_element_pt->bulk_coordinate_derivatives_fct_pt() =
990
991 for (unsigned i = 0; i < 2; i++)
992 {
993 bulk_number = i;
994 face_element_pt->node_pt(i) = node_pt(bulk_number);
995 face_element_pt->bulk_node_number(i) = bulk_number;
996 face_element_pt->normal_sign() = 1;
997 // Set the number of values originally stored at this node
998 face_element_pt->nbulk_value(i) = required_nvalue(bulk_number);
999 // Set the position type for the slope, which is in the s[0]
1000 // direction, so is position_type 1 in "bulk" element
1001 face_element_pt->bulk_position_type(1) = 1;
1002 }
1003 break;
1004 // East face, normal sign is negative
1005 case (1):
1006 // Set the pointer to the bulk coordinate translation scheme
1007 face_element_pt->face_to_bulk_coordinate_fct_pt() =
1009
1010 // Set the pointer to the derivative mappings
1011 face_element_pt->bulk_coordinate_derivatives_fct_pt() =
1013
1014 for (unsigned i = 0; i < 2; i++)
1015 {
1016 bulk_number = 2 * i + 1;
1017 face_element_pt->node_pt(i) = node_pt(bulk_number);
1018 face_element_pt->bulk_node_number(i) = bulk_number;
1019 face_element_pt->normal_sign() = -1;
1020 // Set the number of values originally stored at this node
1021 face_element_pt->nbulk_value(i) = required_nvalue(bulk_number);
1022 // Set the position type for the slope, which is in the s[1]
1023 // direction, so is position_type 2 in the bulk element
1024 face_element_pt->bulk_position_type(1) = 2;
1025 }
1026 break;
1027 // North face, normal sign is negative
1028 case (2):
1029 // Set the pointer to the bulk coordinate translation scheme
1030 face_element_pt->face_to_bulk_coordinate_fct_pt() =
1032
1033 // Set the pointer to the derivative mappings
1034 face_element_pt->bulk_coordinate_derivatives_fct_pt() =
1036
1037 for (unsigned i = 0; i < 2; i++)
1038 {
1039 bulk_number = 2 + i;
1040 face_element_pt->node_pt(i) = node_pt(bulk_number);
1041 face_element_pt->bulk_node_number(i) = bulk_number;
1042 face_element_pt->normal_sign() = -1;
1043 // Set the number of values originally stored at this node
1044 face_element_pt->nbulk_value(i) = required_nvalue(bulk_number);
1045 // Set the position type for the slope, which is in the s[0]
1046 // direction, so is position_type 1 in the "bulk" element
1047 face_element_pt->bulk_position_type(1) = 1;
1048 }
1049 break;
1050
1051 // Now cover the other cases
1052 default:
1053 std::ostringstream error_message;
1054 error_message
1055 << "Face index should only take the values +/- 1 or +/- 2,"
1056 << " not " << face_index << std::endl;
1057 throw OomphLibError(error_message.str(),
1058 OOMPH_CURRENT_FUNCTION,
1059 OOMPH_EXCEPTION_LOCATION);
1060 }
1061 }
1062
1063 /// //////////////////////////////////////////////////////////////////
1064 // 1D SolidQHermiteElements
1065 /// ///////////////////////////////////////////////////////////////////
1066
1067
1068 //=====================================================================
1069 /// Overload the output function
1070 //====================================================================
1071 template<unsigned DIM>
1072 void SolidQHermiteElement<DIM>::output(std::ostream& outfile)
1073 {
1075 }
1076
1077
1078 //=======================================================================
1079 /// The output function for n_plot points in each coordinate direction
1080 /// for the 1D element
1081 //=======================================================================
1082 template<>
1083 void SolidQHermiteElement<1>::output(std::ostream& outfile,
1084 const unsigned& n_plot)
1085 {
1086 // Local variables
1087 Vector<double> s(1);
1088
1089 // Tecplot header info
1090 outfile << "ZONE I=" << n_plot << std::endl;
1091
1092 // Find the dimension of the nodes
1093 unsigned n_dim = this->nodal_dimension();
1094
1095 // Find the Lagrangian dimension of the first node
1096 unsigned n_lagr = static_cast<SolidNode*>(node_pt(0))->nlagrangian();
1097
1098 // Loop over plot points
1099 for (unsigned l = 0; l < n_plot; l++)
1100 {
1101 s[0] = -1.0 + l * 2.0 / (n_plot - 1);
1102
1103 // Output the Eulerian coordinates
1104 for (unsigned i = 0; i < n_dim; i++)
1105 {
1106 outfile << interpolated_x(s, i) << " ";
1107 }
1108
1109 // Output the Lagrangian coordinates
1110 for (unsigned i = 0; i < n_lagr; i++)
1111 {
1112 outfile << interpolated_xi(s, i) << " ";
1113 }
1114 outfile << std::endl;
1115 }
1116 }
1117
1118
1119 //=====================================================================
1120 /// Overload the C-style output function
1121 //====================================================================
1122 template<unsigned DIM>
1124 {
1126 }
1127
1128
1129 //=======================================================================
1130 /// The C-style output function for n_plot points in each coordinate direction
1131 /// for the 1D element
1132 //=======================================================================
1133 template<>
1134 void SolidQHermiteElement<1>::output(FILE* file_pt, const unsigned& n_plot)
1135 {
1136 // Local variables
1137 Vector<double> s(1);
1138
1139 // Tecplot header info
1140 fprintf(file_pt, "ZONE I=%i\n", n_plot);
1141
1142 // Find the dimension of the nodes
1143 unsigned n_dim = this->nodal_dimension();
1144
1145 // Find the Lagrangian dimension of the first node
1146 unsigned n_lagr = static_cast<SolidNode*>(node_pt(0))->nlagrangian();
1147
1148 // Loop over plot points
1149 for (unsigned l = 0; l < n_plot; l++)
1150 {
1151 s[0] = -1.0 + l * 2.0 / (n_plot - 1);
1152
1153 // Output the Eulerian coordinates
1154 for (unsigned i = 0; i < n_dim; i++)
1155 {
1156 fprintf(file_pt, "%g ", interpolated_x(s, i));
1157 }
1158
1159 // Output the Lagrangian coordinates
1160 for (unsigned i = 0; i < n_lagr; i++)
1161 {
1162 fprintf(file_pt, "%g ", interpolated_xi(s, i));
1163 }
1164 fprintf(file_pt, "\n");
1165 }
1166 }
1167
1168
1169 /// ////////////////////////////////////////////////////////////////////////
1170 // 2D SolidQHermiteElements
1171 /// ///////////////////////////////////////////////////////////////////////
1172
1173
1174 //=====================================================================
1175 /// The output function for any number of points per element
1176 //=====================================================================
1177 template<>
1178 void SolidQHermiteElement<2>::output(std::ostream& outfile,
1179 const unsigned& n_p)
1180 {
1181 // Local variables
1182 Vector<double> s(2);
1183
1184 // Tecplot header info
1185 outfile << "ZONE I=" << n_p << ", J=" << n_p << std::endl;
1186
1187 // Find the dimension of the nodes
1188 unsigned n_dim = this->nodal_dimension();
1189
1190 // Find the Lagrangian dimension of the first node
1191 unsigned n_lagr = static_cast<SolidNode*>(node_pt(0))->nlagrangian();
1192
1193 // Loop over element nodes
1194 for (unsigned l2 = 0; l2 < n_p; l2++)
1195 {
1196 s[1] = -1.0 + l2 * 2.0 / (n_p - 1);
1197 for (unsigned l1 = 0; l1 < n_p; l1++)
1198 {
1199 s[0] = -1.0 + l1 * 2.0 / (n_p - 1);
1200
1201 // Output the Eulerian coordinates
1202 for (unsigned i = 0; i < n_dim; i++)
1203 {
1204 outfile << interpolated_x(s, i) << " ";
1205 }
1206
1207 // Output the Lagrangian coordinates
1208 for (unsigned i = 0; i < n_lagr; i++)
1209 {
1210 outfile << interpolated_xi(s, i) << " ";
1211 }
1212 outfile << std::endl;
1213 }
1214 }
1215 outfile << std::endl;
1216 }
1217
1218
1219 //=====================================================================
1220 /// The C-style output function for any number of points per element
1221 //=====================================================================
1222 template<>
1223 void SolidQHermiteElement<2>::output(FILE* file_pt, const unsigned& n_plot)
1224 {
1225 // Local variables
1226 Vector<double> s(2);
1227
1228 // Tecplot header info
1229 fprintf(file_pt, "ZONE I=%i, J=%i\n", n_plot, n_plot);
1230
1231 // Find the dimension of the nodes
1232 unsigned n_dim = this->nodal_dimension();
1233
1234 // Find the Lagrangian dimension of the first node
1235 unsigned n_lagr = static_cast<SolidNode*>(node_pt(0))->nlagrangian();
1236
1237 // Loop over element nodes
1238 for (unsigned l2 = 0; l2 < n_plot; l2++)
1239 {
1240 s[1] = -1.0 + l2 * 2.0 / (n_plot - 1);
1241 for (unsigned l1 = 0; l1 < n_plot; l1++)
1242 {
1243 s[0] = -1.0 + l1 * 2.0 / (n_plot - 1);
1244
1245 // Output the Eulerian coordinates
1246 for (unsigned i = 0; i < n_dim; i++)
1247 {
1248 fprintf(file_pt, "%g ", interpolated_x(s, i));
1249 }
1250
1251 // Output the Lagrangian coordinates
1252 for (unsigned i = 0; i < n_lagr; i++)
1253 {
1254 fprintf(file_pt, "%g ", interpolated_xi(s, i));
1255 }
1256 fprintf(file_pt, "\n");
1257 }
1258 }
1259 fprintf(file_pt, "\n");
1260 }
1261
1262
1263 //=======================================================================
1264 /// Function to setup geometrical information for lower-dimensional
1265 /// FaceElements for the solid hermite elements. We need to
1266 /// construct the basic element and then sort out the Lagrangian
1267 /// coordinates
1268 //=======================================================================
1269 template<unsigned DIM>
1271 const int& face_index, FaceElement* face_element_pt)
1272 {
1273 // Build the standard non-solid FaceElement
1274 QHermiteElement<DIM>::build_face_element(face_index, face_element_pt);
1275
1276 // Set the Lagrangian dimension from the first node of the present element
1277 dynamic_cast<SolidFiniteElement*>(face_element_pt)
1278 ->set_lagrangian_dimension(
1279 static_cast<SolidNode*>(node_pt(0))->nlagrangian());
1280 }
1281
1282 //==================================================================
1283 /// Force build of templates
1284 //==================================================================
1285 template class QHermiteElement<1>;
1286 template class QHermiteElement<2>;
1287 template class DiagQHermiteElement<1>;
1288 template class DiagQHermiteElement<2>;
1289
1290
1291 template class SolidQHermiteElement<1>;
1292 template class SolidQHermiteElement<2>;
1293 template class SolidDiagQHermiteElement<1>;
1294 template class SolidDiagQHermiteElement<2>;
1295
1296} // namespace oomph
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
Definition: shape.h:278
These elements are exactly the same as QHermiteElements, but they employ the simplifying assumption t...
FaceElements are elements that coincide with the faces of higher-dimensional "bulk" elements....
Definition: elements.h:4338
int & normal_sign()
Sign of outer unit normal (relative to cross-products of tangent vectors in the corresponding "bulk" ...
Definition: elements.h:4612
int & face_index()
Index of the face (a number that uniquely identifies the face in the element)
Definition: elements.h:4626
void bulk_node_number_resize(const unsigned &i)
Resize the storage for the bulk node numbers.
Definition: elements.h:4818
BulkCoordinateDerivativesFctPt & bulk_coordinate_derivatives_fct_pt()
Return the pointer to the function that returns the derivatives of the bulk coordinates wrt the face ...
Definition: elements.h:4770
unsigned & nbulk_value(const unsigned &n)
Return the number of values originally stored at local node n (before the FaceElement added additiona...
Definition: elements.h:4845
unsigned & bulk_node_number(const unsigned &n)
Return the bulk node number that corresponds to the n-th local node number.
Definition: elements.h:4825
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
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
CoordinateMappingFctPt & face_to_bulk_coordinate_fct_pt()
Return the pointer to the function that maps the face coordinate to the bulk coordinate.
Definition: elements.h:4755
void nbulk_value_resize(const unsigned &i)
Resize the storage for the number of values originally stored at the local nodes to i entries.
Definition: elements.h:4860
void set_nodal_dimension(const unsigned &nodal_dim)
Set the dimension of the nodes in the element. This will typically only be required when constructing...
Definition: elements.h:1390
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
void set_halo(const unsigned &non_halo_proc_ID)
Label the element as halo and specify processor that holds non-halo counterpart.
Definition: elements.h:1151
An OomphLibError object which should be thrown when an run-time error is encountered....
/////////////////////////////////////////////////////////////////// /////////////////////////////////...
void shape(const Vector< double > &s, Shape &psi) const
Function to calculate the geometric shape functions at local coordinate s.
static Gauss< DIM, 3 > Default_integration_scheme
Default integration rule: Gaussian integration of same 'order' as the element.
void build_face_element(const int &face_index, FaceElement *face_element_pt)
Build the lower-dimensional FaceElement of the type QHermiteElement<DIM-1>. The face index takes a va...
void dshape_local(const Vector< double > &s, Shape &psi, DShape &dpsids) const
Function to compute the geometric shape functions and derivatives w.r.t. local coordinates at local c...
void output(std::ostream &outfile)
Output.
void d2shape_local(const Vector< double > &s, Shape &psi, DShape &dpsids, DShape &d2psids) const
Function to compute the geometric shape functions and also first and second derivatives wrt local coo...
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition: elements.h:3561
A Class for nodes that deform elastically (i.e. position is an unknown in the problem)....
Definition: nodes.h:1686
unsigned nlagrangian() const
Return number of lagrangian coordinates.
Definition: nodes.h:1870
////////////////////////////////////////////////////////////////////
void output(std::ostream &outfile)
Overload the output function.
void output()
Doc the command line arguments.
void dshape(const double &s, double DPsi[2][2])
Derivatives of 1D Hermite shape functions.
Definition: shape.h:1206
void d2shape(const double &s, double DPsi[2][2])
Second derivatives of the Hermite shape functions.
Definition: shape.h:1217
void shape(const double &s, double Psi[2][2])
Constructor sets the values of the shape functions at the position s.
Definition: shape.h:1194
void faces0(const Vector< double > &s, DenseMatrix< double > &dsbulk_dsface, unsigned &interior_direction)
Function for both faces – the bulk coordinate is fixed on both.
void face1(const Vector< double > &s, Vector< double > &s_bulk)
The translation scheme for the face s0 = 1.0.
void face0(const Vector< double > &s, Vector< double > &s_bulk)
The translation scheme for the face s0 = -1.0.
void faces1(const Vector< double > &s, DenseMatrix< double > &dsbulk_dsface, unsigned &interior_direction)
Function for the north and south faces, along which s1 is fixed.
void faces0(const Vector< double > &s, DenseMatrix< double > &dsbulk_dsface, unsigned &interior_direction)
Function for the east and west faces, along which s0 is fixed.
void face2(const Vector< double > &s, Vector< double > &s_bulk)
The translation scheme for the east face (s0 = 1.0)
void face1(const Vector< double > &s, Vector< double > &s_bulk)
The translation scheme for the south face (s1 = -1.0)
void face0(const Vector< double > &s, Vector< double > &s_bulk)
The translation scheme for the west face (s0 = -1.0)
void face3(const Vector< double > &s, Vector< double > &s_bulk)
The translation scheme for the north face (s1 = 1.0)
//////////////////////////////////////////////////////////////////// ////////////////////////////////...