Qelements.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 member functions for Qelements
27
28// Include the appropriate headers
29#include "Qelements.h"
30
31namespace oomph
32{
33 /// /////////////////////////////////////////////////////////////
34 // 1D Qelements
35 /// /////////////////////////////////////////////////////////////
36
37 //=======================================================================
38 /// Assign the static Default_integration_scheme
39 //=======================================================================
40 template<unsigned NNODE_1D>
42
43
44 //==================================================================
45 /// Return the node at the specified local coordinate
46 //==================================================================
47 template<unsigned NNODE_1D>
49 const Vector<double>& s) const
50 {
51 // Load the tolerance into a local variable
53 // There is one possible index.
54 Vector<int> index(1);
55
56 // Determine the index
57 // -------------------
58
59 // If we are at the lower limit, the index is zero
60 if (std::fabs(s[0] + 1.0) < tol)
61 {
62 index[0] = 0;
63 }
64 // If we are at the upper limit, the index is the number of nodes minus 1
65 else if (std::fabs(s[0] - 1.0) < tol)
66 {
67 index[0] = NNODE_1D - 1;
68 }
69 // Otherwise, we have to calculate the index in general
70 else
71 {
72 // For uniformly spaced nodes the node number would be
73 double float_index = 0.5 * (1.0 + s[0]) * (NNODE_1D - 1);
74 // Convert to an integer by taking the floor (rounding down)
75 index[0] = static_cast<int>(std::floor(float_index));
76 // What is the excess. This should be safe because the
77 // we have rounded down
78 double excess = float_index - index[0];
79 // If the excess is bigger than our tolerance there is no node,
80 // return null
81 // Note that we test at both lower and upper ends.
82 if ((excess > tol) && ((1.0 - excess) > tol))
83 {
84 return 0;
85 }
86 // If we are at the upper end (i.e. the system has actually rounded up)
87 // we need to add one to the index
88 if ((1.0 - excess) <= tol)
89 {
90 index[0] += 1;
91 }
92 }
93 // If we've got here we have a node, so let's return a pointer to it
94 return node_pt(index[0]);
95 }
96
97
98 //=======================================================================
99 /// Shape function for specific QElement<1,..>
100 //=======================================================================
101 template<unsigned NNODE_1D>
103 {
104 // Local storage for the shape functions
105 double Psi[NNODE_1D];
106 // Call the OneDimensional Shape functions
107 OneDimLagrange::shape<NNODE_1D>(s[0], Psi);
108
109 // Now let's loop over the nodal points in the element
110 // and copy the values back in
111 for (unsigned i = 0; i < NNODE_1D; i++)
112 {
113 psi[i] = Psi[i];
114 }
115 }
116
117 //=======================================================================
118 /// Derivatives of shape functions for specific QElement<1,..>
119 //=======================================================================
120 template<unsigned NNODE_1D>
122 Shape& psi,
123 DShape& dpsids) const
124 {
125 // Local storage
126 double Psi[NNODE_1D];
127 double DPsi[NNODE_1D];
128 // Call the shape functions and derivatives
129 OneDimLagrange::shape<NNODE_1D>(s[0], Psi);
130 OneDimLagrange::dshape<NNODE_1D>(s[0], DPsi);
131
132 // Loop over shape functions in element and assign to psi
133 for (unsigned l = 0; l < NNODE_1D; l++)
134 {
135 psi[l] = Psi[l];
136 dpsids(l, 0) = DPsi[l];
137 }
138 }
139
140 //=======================================================================
141 /// Second derivatives of shape functions for specific QElement<1,..>:
142 /// d2psids(i,0) = \f$ d^2 \psi_j / d s^2 \f$
143 //=======================================================================
144 template<unsigned NNODE_1D>
146 Shape& psi,
147 DShape& dpsids,
148 DShape& d2psids) const
149 {
150 // Local storage for the shape functions
151 double Psi[NNODE_1D];
152 double DPsi[NNODE_1D];
153 double D2Psi[NNODE_1D];
154 // Call the shape functions and derivatives
155 OneDimLagrange::shape<NNODE_1D>(s[0], Psi);
156 OneDimLagrange::dshape<NNODE_1D>(s[0], DPsi);
157 OneDimLagrange::d2shape<NNODE_1D>(s[0], D2Psi);
158
159 // Loop over shape functions in element and assign to psi
160 for (unsigned l = 0; l < NNODE_1D; l++)
161 {
162 psi[l] = Psi[l];
163 dpsids(l, 0) = DPsi[l];
164 d2psids(l, 0) = D2Psi[l];
165 }
166 }
167
168 //=======================================================================
169 /// The output function for general 1D QElements
170 //=======================================================================
171 template<unsigned NNODE_1D>
172 void QElement<1, NNODE_1D>::output(std::ostream& outfile)
173 {
174 output(outfile, NNODE_1D);
175 }
176
177 //=======================================================================
178 /// The output function for n_plot points in each coordinate direction
179 //=======================================================================
180 template<unsigned NNODE_1D>
181 void QElement<1, NNODE_1D>::output(std::ostream& outfile,
182 const unsigned& n_plot)
183 {
184 // Local variables
185 Vector<double> s(1);
186
187 // Tecplot header info
188 outfile << "ZONE I=" << n_plot << std::endl;
189
190 // Find the dimension of the nodes
191 unsigned n_dim = this->nodal_dimension();
192
193 // Loop over plot points
194 for (unsigned l = 0; l < n_plot; l++)
195 {
196 s[0] = -1.0 + l * 2.0 / (n_plot - 1);
197 // Output the coordinates
198 for (unsigned i = 0; i < n_dim; i++)
199 {
200 outfile << interpolated_x(s, i) << " ";
201 }
202 outfile << std::endl;
203 }
204 outfile << std::endl;
205 }
206
207
208 //=======================================================================
209 /// C style output function for general 1D QElements
210 //=======================================================================
211 template<unsigned NNODE_1D>
213 {
214 output(file_pt, NNODE_1D);
215 }
216
217 //=======================================================================
218 /// C style output function for n_plot points in each coordinate direction
219 //=======================================================================
220 template<unsigned NNODE_1D>
221 void QElement<1, NNODE_1D>::output(FILE* file_pt, const unsigned& n_plot)
222 {
223 // Local variables
224 Vector<double> s(1);
225
226 // Tecplot header info
227 // outfile << "ZONE I=" << n_plot << std::endl;
228 fprintf(file_pt, "ZONE I=%i\n", n_plot);
229
230 // Find the dimension of the first node
231 unsigned n_dim = this->nodal_dimension();
232
233 // Loop over plot points
234 for (unsigned l = 0; l < n_plot; l++)
235 {
236 s[0] = -1.0 + l * 2.0 / (n_plot - 1);
237
238 // Output the coordinates
239 for (unsigned i = 0; i < n_dim; i++)
240 {
241 // outfile << interpolated_x(s,i) << " " ;
242 fprintf(file_pt, "%g ", interpolated_x(s, i));
243 }
244 // outfile << std::endl;
245 fprintf(file_pt, "\n");
246 }
247 // outfile << std::endl;
248 fprintf(file_pt, "\n");
249 }
250
251
252 /// /////////////////////////////////////////////////////////////
253 // 2D Qelements
254 /// /////////////////////////////////////////////////////////////
255
256 /// Assign the spatial integration scheme
257 template<unsigned NNODE_1D>
259
260
261 //==================================================================
262 /// Return the node at the specified local coordinate
263 //==================================================================
264 template<unsigned NNODE_1D>
266 const Vector<double>& s) const
267 {
268 // Load the tolerance into a local variable
270 // There are two possible indices.
271 Vector<int> index(2);
272 // Loop over the coordinates and determine the indices
273 for (unsigned i = 0; i < 2; i++)
274 {
275 // If we are at the lower limit, the index is zero
276 if (std::fabs(s[i] + 1.0) < tol)
277 {
278 index[i] = 0;
279 }
280 // If we are at the upper limit, the index is the number of nodes minus 1
281 else if (std::fabs(s[i] - 1.0) < tol)
282 {
283 index[i] = NNODE_1D - 1;
284 }
285 // Otherwise, we have to calculate the index in general
286 else
287 {
288 // For uniformly spaced nodes the node number would be
289 double float_index = 0.5 * (1.0 + s[i]) * (NNODE_1D - 1);
290 // Convert to an integer by taking the floor (rounding down)
291 index[i] = static_cast<int>(std::floor(float_index));
292 // What is the excess. This should be safe because the
293 // we have rounded down
294 double excess = float_index - index[i];
295 // If the excess is bigger than our tolerance there is no node,
296 // return null
297 // Note that we test at both lower and upper ends.
298 if ((excess > tol) && ((1.0 - excess) > tol))
299 {
300 return 0;
301 }
302 // If we are at the upper end (i.e. the system has actually rounded up)
303 // we need to add one to the index
304 if ((1.0 - excess) <= tol)
305 {
306 index[i] += 1;
307 }
308 }
309 }
310 // If we've got here we have a node, so let's return a pointer to it
311 return node_pt(index[0] + NNODE_1D * index[1]);
312 }
313
314
315 //=======================================================================
316 /// Shape function for specific QElement<2,..>
317 ///
318 //=======================================================================
319 template<unsigned NNODE_1D>
321 {
322 // Local storage
323 double Psi[2][NNODE_1D];
324 // Call the OneDimensional Shape functions
325 OneDimLagrange::shape<NNODE_1D>(s[0], Psi[0]);
326 OneDimLagrange::shape<NNODE_1D>(s[1], Psi[1]);
327 // Index for the shape functions
328 unsigned index = 0;
329
330 // Now let's loop over the nodal points in the element
331 // s1 is the "x" coordinate, s2 the "y"
332 for (unsigned i = 0; i < NNODE_1D; i++)
333 {
334 for (unsigned j = 0; j < NNODE_1D; j++)
335 {
336 // Multiply the two 1D functions together to get the 2D function
337 psi[index] = Psi[1][i] * Psi[0][j];
338 // Incremenet the index
339 ++index;
340 }
341 }
342 }
343
344 //=======================================================================
345 /// Derivatives of shape functions for specific QElement<2,..>
346 //=======================================================================
347 template<unsigned NNODE_1D>
349 Shape& psi,
350 DShape& dpsids) const
351 {
352 // Local storage
353 double Psi[2][NNODE_1D];
354 double DPsi[2][NNODE_1D];
355 unsigned index = 0;
356
357 // Call the shape functions and derivatives
358 OneDimLagrange::shape<NNODE_1D>(s[0], Psi[0]);
359 OneDimLagrange::shape<NNODE_1D>(s[1], Psi[1]);
360 OneDimLagrange::dshape<NNODE_1D>(s[0], DPsi[0]);
361 OneDimLagrange::dshape<NNODE_1D>(s[1], DPsi[1]);
362
363 // Loop over shape functions in element
364 for (unsigned i = 0; i < NNODE_1D; i++)
365 {
366 for (unsigned j = 0; j < NNODE_1D; j++)
367 {
368 // Assign the values
369 dpsids(index, 0) = Psi[1][i] * DPsi[0][j];
370 dpsids(index, 1) = DPsi[1][i] * Psi[0][j];
371 psi[index] = Psi[1][i] * Psi[0][j];
372 // Increment the index
373 ++index;
374 }
375 }
376 }
377
378
379 //=======================================================================
380 /// Second derivatives of shape functions for specific QElement<2,..>:
381 /// d2psids(i,0) = \f$ \partial^2 \psi_j / \partial s_0^2 \f$
382 /// d2psids(i,1) = \f$ \partial^2 \psi_j / \partial s_1^2 \f$
383 /// d2psids(i,2) = \f$ \partial^2 \psi_j / \partial s_0 \partial s_1 \f$
384 //=======================================================================
385 template<unsigned NNODE_1D>
387 Shape& psi,
388 DShape& dpsids,
389 DShape& d2psids) const
390 {
391 // Local storage
392 double Psi[2][NNODE_1D];
393 double DPsi[2][NNODE_1D];
394 double D2Psi[2][NNODE_1D];
395 // Index for the assembly process
396 unsigned index = 0;
397
398 // Call the shape functions and derivatives
399 OneDimLagrange::shape<NNODE_1D>(s[0], Psi[0]);
400 OneDimLagrange::shape<NNODE_1D>(s[1], Psi[1]);
401 OneDimLagrange::dshape<NNODE_1D>(s[0], DPsi[0]);
402 OneDimLagrange::dshape<NNODE_1D>(s[1], DPsi[1]);
403 OneDimLagrange::d2shape<NNODE_1D>(s[0], D2Psi[0]);
404 OneDimLagrange::d2shape<NNODE_1D>(s[1], D2Psi[1]);
405
406 // Loop over shape functions in element
407 for (unsigned i = 0; i < NNODE_1D; i++)
408 {
409 for (unsigned j = 0; j < NNODE_1D; j++)
410 {
411 // Assign the values
412 psi[index] = Psi[1][i] * Psi[0][j];
413 // First derivatives
414 dpsids(index, 0) = Psi[1][i] * DPsi[0][j];
415 dpsids(index, 1) = DPsi[1][i] * Psi[0][j];
416 // Second derivatives
417 // N.B. index 2 is the mixed derivative
418 d2psids(index, 0) = Psi[1][i] * D2Psi[0][j];
419 d2psids(index, 1) = D2Psi[1][i] * Psi[0][j];
420 d2psids(index, 2) = DPsi[1][i] * DPsi[0][j];
421 // Increment the index
422 ++index;
423 }
424 }
425 }
426
427
428 //===========================================================
429 /// The output function for QElement<2,NNODE_1D>
430 //===========================================================
431 template<unsigned NNODE_1D>
432 void QElement<2, NNODE_1D>::output(std::ostream& outfile)
433 {
434 output(outfile, NNODE_1D);
435 }
436
437 //=======================================================================
438 /// The output function for n_plot points in each coordinate direction
439 //=======================================================================
440 template<unsigned NNODE_1D>
441 void QElement<2, NNODE_1D>::output(std::ostream& outfile,
442 const unsigned& n_plot)
443 {
444 // Local variables
445 Vector<double> s(2);
446
447 // Tecplot header info
448 outfile << "ZONE I=" << n_plot << ", J=" << n_plot << std::endl;
449
450 // Find the dimension of the nodes
451 unsigned n_dim = this->nodal_dimension();
452
453 // Loop over plot points
454 for (unsigned l2 = 0; l2 < n_plot; l2++)
455 {
456 s[1] = -1.0 + l2 * 2.0 / (n_plot - 1);
457 for (unsigned l1 = 0; l1 < n_plot; l1++)
458 {
459 s[0] = -1.0 + l1 * 2.0 / (n_plot - 1);
460
461 // Output the coordinates
462 for (unsigned i = 0; i < n_dim; i++)
463 {
464 outfile << interpolated_x(s, i) << " ";
465 }
466 outfile << std::endl;
467 }
468 }
469 outfile << std::endl;
470 }
471
472
473 //===========================================================
474 /// C-style output function for QElement<2,NNODE_1D>
475 //===========================================================
476 template<unsigned NNODE_1D>
478 {
479 output(file_pt, NNODE_1D);
480 }
481
482 //=======================================================================
483 /// C-style output function for n_plot points in each coordinate direction
484 //=======================================================================
485 template<unsigned NNODE_1D>
486 void QElement<2, NNODE_1D>::output(FILE* file_pt, const unsigned& n_plot)
487 {
488 // Local variables
489 Vector<double> s(2);
490
491 // Find the dimension of the nodes
492 unsigned n_dim = this->nodal_dimension();
493
494 // Tecplot header info
495 // outfile << "ZONE I=" << n_plot << ", J=" << n_plot << std::endl;
496 fprintf(file_pt, "ZONE I=%i, J=%i\n", n_plot, n_plot);
497
498 // Loop over element nodes
499 for (unsigned l2 = 0; l2 < n_plot; l2++)
500 {
501 s[1] = -1.0 + l2 * 2.0 / (n_plot - 1);
502 for (unsigned l1 = 0; l1 < n_plot; l1++)
503 {
504 s[0] = -1.0 + l1 * 2.0 / (n_plot - 1);
505
506 // Output the coordinates
507 for (unsigned i = 0; i < n_dim; i++)
508 {
509 // outfile << interpolated_x(s,i) << " " ;
510 fprintf(file_pt, "%g ", interpolated_x(s, i));
511 }
512 // outfile << std::endl;
513 fprintf(file_pt, "\n");
514 }
515 }
516 // outfile << std::endl;
517 fprintf(file_pt, "\n");
518 }
519
520
521 /// /////////////////////////////////////////////////////////////
522 // 3D Qelements
523 /// /////////////////////////////////////////////////////////////
524
525 /// Assign the spatial integration scheme
526 template<unsigned NNODE_1D>
528
529 //==================================================================
530 /// Return the node at the specified local coordinate
531 //==================================================================
532 template<unsigned NNODE_1D>
534 const Vector<double>& s) const
535 {
536 // Load the tolerance into a local variable
538 // There are now three possible indices
539 Vector<int> index(3);
540 // Loop over the coordinates
541 for (unsigned i = 0; i < 3; i++)
542 {
543 // If we are at the lower limit, the index is zero
544 if (std::fabs(s[i] + 1.0) < tol)
545 {
546 index[i] = 0;
547 }
548 // If we are at the upper limit, the index is the number of nodes minus 1
549 else if (std::fabs(s[i] - 1.0) < tol)
550 {
551 index[i] = NNODE_1D - 1;
552 }
553 // Otherwise, we have to calculate the index in general
554 else
555 {
556 // For uniformly spaced nodes the node number would be
557 double float_index = 0.5 * (1.0 + s[i]) * (NNODE_1D - 1);
558 // Conver to an integer by taking the floor (rounding down)
559 index[i] = static_cast<int>(std::floor(float_index));
560 // What is the excess. This should be safe because
561 // we have rounded down
562 double excess = float_index - index[i];
563 // If the excess is bigger than our tolerance there is no node,
564 // return null. Note that we test at both ends
565 if ((excess > tol) && ((1.0 - excess) > tol))
566 {
567 return 0;
568 }
569 // If we are at the upper end (i.e. the system has actually rounded up)
570 // we need to add one to the index
571 if ((1.0 - excess) <= tol)
572 {
573 index[i] += 1;
574 }
575 }
576 }
577 // If we've got here we have a node, so let's return a pointer to it
578 return node_pt(index[0] + NNODE_1D * index[1] +
579 NNODE_1D * NNODE_1D * index[2]);
580 }
581
582
583 //=======================================================================
584 /// Shape function for specific QElement<3,..>
585 //=======================================================================
586 template<unsigned NNODE_1D>
588 {
589 // Local storage
590 double Psi[3][NNODE_1D];
591
592 // Call the OneDimensional Shape functions
593 OneDimLagrange::shape<NNODE_1D>(s[0], Psi[0]);
594 OneDimLagrange::shape<NNODE_1D>(s[1], Psi[1]);
595 OneDimLagrange::shape<NNODE_1D>(s[2], Psi[2]);
596
597 // Index for the shape functions
598 unsigned index = 0;
599
600 // Now let's loop over the nodal points in the element
601 // s1 is the "x" coordinate, s2 the "y"
602 for (unsigned i = 0; i < NNODE_1D; i++)
603 {
604 for (unsigned j = 0; j < NNODE_1D; j++)
605 {
606 for (unsigned k = 0; k < NNODE_1D; k++)
607 {
608 /*Multiply the three 1D functions together to get the 3D function*/
609 psi[index] = Psi[2][i] * Psi[1][j] * Psi[0][k];
610 // Increment the index
611 ++index;
612 }
613 }
614 }
615 }
616
617 //=======================================================================
618 /// Derivatives of shape functions for specific QElement<3,..>
619 //=======================================================================
620 template<unsigned NNODE_1D>
622 Shape& psi,
623 DShape& dpsids) const
624 {
625 // Local storage
626 double Psi[3][NNODE_1D];
627 double DPsi[3][NNODE_1D];
628 // Index of the total shape function
629 unsigned index = 0;
630
631 // Call the OneDimensional Shape functions and derivatives
632 OneDimLagrange::shape<NNODE_1D>(s[0], Psi[0]);
633 OneDimLagrange::shape<NNODE_1D>(s[1], Psi[1]);
634 OneDimLagrange::shape<NNODE_1D>(s[2], Psi[2]);
635 OneDimLagrange::dshape<NNODE_1D>(s[0], DPsi[0]);
636 OneDimLagrange::dshape<NNODE_1D>(s[1], DPsi[1]);
637 OneDimLagrange::dshape<NNODE_1D>(s[2], DPsi[2]);
638
639
640 // Loop over shape functions in element
641 for (unsigned i = 0; i < NNODE_1D; i++)
642 {
643 for (unsigned j = 0; j < NNODE_1D; j++)
644 {
645 for (unsigned k = 0; k < NNODE_1D; k++)
646 {
647 // Assign the values
648 dpsids(index, 0) = Psi[2][i] * Psi[1][j] * DPsi[0][k];
649 dpsids(index, 1) = Psi[2][i] * DPsi[1][j] * Psi[0][k];
650 dpsids(index, 2) = DPsi[2][i] * Psi[1][j] * Psi[0][k];
651
652 psi[index] = Psi[2][i] * Psi[1][j] * Psi[0][k];
653 // Increment the index
654 ++index;
655 }
656 }
657 }
658 }
659
660 //=======================================================================
661 /// Second derivatives of shape functions for specific QElement<3,..>:
662 /// d2psids(i,0) = \f$ \partial^2 \psi_j / \partial s_0^2 \f$
663 /// d2psids(i,1) = \f$ \partial^2 \psi_j / \partial s_1^2 \f$
664 /// d2psids(i,2) = \f$ \partial^2 \psi_j / \partial s_2^2 \f$
665 /// d2psids(i,3) = \f$ \partial^2 \psi_j / \partial s_0 \partial s_1 \f$
666 /// d2psids(i,4) = \f$ \partial^2 \psi_j / \partial s_0 \partial s_2 \f$
667 /// d2psids(i,5) = \f$ \partial^2 \psi_j / \partial s_1 \partial s_2 \f$
668 //=======================================================================
669 template<unsigned NNODE_1D>
671 Shape& psi,
672 DShape& dpsids,
673 DShape& d2psids) const
674 {
675 // Local storage
676 double Psi[3][NNODE_1D];
677 double DPsi[3][NNODE_1D];
678 double D2Psi[3][NNODE_1D];
679 // Index of the shape function
680 unsigned index = 0;
681
682 // Call the OneDimensional Shape functions and derivatives
683 OneDimLagrange::shape<NNODE_1D>(s[0], Psi[0]);
684 OneDimLagrange::shape<NNODE_1D>(s[1], Psi[1]);
685 OneDimLagrange::shape<NNODE_1D>(s[2], Psi[2]);
686 OneDimLagrange::dshape<NNODE_1D>(s[0], DPsi[0]);
687 OneDimLagrange::dshape<NNODE_1D>(s[1], DPsi[1]);
688 OneDimLagrange::dshape<NNODE_1D>(s[2], DPsi[2]);
689 OneDimLagrange::d2shape<NNODE_1D>(s[0], D2Psi[0]);
690 OneDimLagrange::d2shape<NNODE_1D>(s[1], D2Psi[1]);
691 OneDimLagrange::d2shape<NNODE_1D>(s[2], D2Psi[2]);
692
693 // Loop over shape functions in element
694 for (unsigned i = 0; i < NNODE_1D; i++)
695 {
696 for (unsigned j = 0; j < NNODE_1D; j++)
697 {
698 for (unsigned k = 0; k < NNODE_1D; k++)
699 {
700 // Assign the values
701 psi[index] = Psi[2][i] * Psi[1][j] * Psi[0][k];
702
703 dpsids(index, 0) = Psi[2][i] * Psi[1][j] * DPsi[0][k];
704 dpsids(index, 1) = Psi[2][i] * DPsi[1][j] * Psi[0][k];
705 dpsids(index, 2) = DPsi[2][i] * Psi[1][j] * Psi[0][k];
706
707 // Second derivative values
708 d2psids(index, 0) = Psi[2][i] * Psi[1][j] * D2Psi[0][k];
709 d2psids(index, 1) = Psi[2][i] * D2Psi[1][j] * Psi[0][k];
710 d2psids(index, 2) = D2Psi[2][i] * Psi[1][j] * Psi[0][k];
711 // Convention for higher indices
712 // 3: mixed 12, 4: mixed 13, 5: mixed 23
713 d2psids(index, 3) = Psi[2][i] * DPsi[1][j] * DPsi[0][k];
714 d2psids(index, 4) = DPsi[2][i] * Psi[1][j] * DPsi[0][k];
715 d2psids(index, 5) = DPsi[2][i] * DPsi[1][j] * Psi[0][k];
716 // Increment the index
717 ++index;
718 }
719 }
720 }
721 }
722
723 //===========================================================
724 /// The output function for QElement<3,NNODE_1D>
725 //===========================================================
726 template<unsigned NNODE_1D>
727 void QElement<3, NNODE_1D>::output(std::ostream& outfile)
728 {
729 output(outfile, NNODE_1D);
730 }
731
732 //=======================================================================
733 /// The output function for n_plot points in each coordinate direction
734 //=======================================================================
735 template<unsigned NNODE_1D>
736 void QElement<3, NNODE_1D>::output(std::ostream& outfile,
737 const unsigned& n_plot)
738 {
739 // Local variables
740 Vector<double> s(3);
741
742 // Tecplot header info
743 outfile << "ZONE I=" << n_plot << ", J=" << n_plot << ", K=" << n_plot
744 << std::endl;
745
746 // Find the dimension of the nodes
747 unsigned n_dim = this->nodal_dimension();
748
749 // Loop over element nodes
750 for (unsigned l3 = 0; l3 < n_plot; l3++)
751 {
752 s[2] = -1.0 + l3 * 2.0 / (n_plot - 1);
753 for (unsigned l2 = 0; l2 < n_plot; l2++)
754 {
755 s[1] = -1.0 + l2 * 2.0 / (n_plot - 1);
756 for (unsigned l1 = 0; l1 < n_plot; l1++)
757 {
758 s[0] = -1.0 + l1 * 2.0 / (n_plot - 1);
759
760 // Output the coordinates
761 for (unsigned i = 0; i < n_dim; i++)
762 {
763 outfile << interpolated_x(s, i) << " ";
764 }
765 outfile << std::endl;
766 }
767 }
768 }
769 outfile << std::endl;
770 }
771
772
773 //===========================================================
774 /// C-style output function for QElement<3,NNODE_1D>
775 //===========================================================
776 template<unsigned NNODE_1D>
778 {
779 output(file_pt, NNODE_1D);
780 }
781
782 //=======================================================================
783 /// C-style output function for n_plot points in each coordinate direction
784 //=======================================================================
785 template<unsigned NNODE_1D>
786 void QElement<3, NNODE_1D>::output(FILE* file_pt, const unsigned& n_plot)
787 {
788 // Local variables
789 Vector<double> s(3);
790
791 // Tecplot header info
792 fprintf(file_pt, "ZONE I=%i, J=%i, K=%i\n", n_plot, n_plot, n_plot);
793
794 // Find the dimension of the nodes
795 unsigned n_dim = this->nodal_dimension();
796
797 // Loop over element nodes
798 for (unsigned l3 = 0; l3 < n_plot; l3++)
799 {
800 s[2] = -1.0 + l3 * 2.0 / (n_plot - 1);
801 for (unsigned l2 = 0; l2 < n_plot; l2++)
802 {
803 s[1] = -1.0 + l2 * 2.0 / (n_plot - 1);
804 for (unsigned l1 = 0; l1 < n_plot; l1++)
805 {
806 s[0] = -1.0 + l1 * 2.0 / (n_plot - 1);
807
808 // Output the coordinates
809 for (unsigned i = 0; i < n_dim; i++)
810 {
811 fprintf(file_pt, "%g ", interpolated_x(s, i));
812 }
813 fprintf(file_pt, "\n");
814 }
815 }
816 }
817 fprintf(file_pt, "\n");
818 }
819
820
821 //===================================================================
822 // Build required templates
823 //===================================================================
824 template class QElement<1, 2>;
825 template class QElement<1, 3>;
826 template class QElement<1, 4>;
827
828 template class QElement<2, 2>;
829 template class QElement<2, 3>;
830 template class QElement<2, 4>;
831
832 template class QElement<3, 2>;
833 template class QElement<3, 3>;
834 template class QElement<3, 4>;
835
836} // 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
static const double Node_location_tolerance
Default value for the tolerance to be used when locating nodes via local coordinates.
Definition: elements.h:1374
Class for multidimensional Gaussian integration rules.
Definition: integral.h:145
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
static Gauss< 1, NNODE_1D > Default_integration_scheme
Default integration rule: Gaussian integration of same 'order' as the element.
Definition: Qelements.h:490
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
Definition: Qelements.h:459
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
void output()
Doc the command line arguments.
void shape(const double &s, double *Psi)
Definition for 1D Lagrange shape functions. The value of all the shape functions at the local coordin...
Definition: shape.h:564
//////////////////////////////////////////////////////////////////// ////////////////////////////////...