Telements.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 Telements
27
28
29// oomph-lib headers
30#include "Telements.h"
31
32
33namespace oomph
34{
35 //=======================================================================
36 /// Assign the static integral
37 //=======================================================================
38 template<unsigned NNODE_1D>
40 template<unsigned NNODE_1D>
42 template<unsigned NNODE_1D>
44
45 /// ///////////////////////////////////////////////////////////////
46 // 1D Telements
47 /// ///////////////////////////////////////////////////////////////
48
49 //=======================================================================
50 /// The output function for general 1D TElements
51 //=======================================================================
52 template<unsigned NNODE_1D>
53 void TElement<1, NNODE_1D>::output(std::ostream& outfile)
54 {
55 output(outfile, NNODE_1D);
56 }
57
58 //=======================================================================
59 /// The output function for n_plot points in each coordinate direction
60 //=======================================================================
61 template<unsigned NNODE_1D>
62 void TElement<1, NNODE_1D>::output(std::ostream& outfile,
63 const unsigned& n_plot)
64 {
65 // Local variables
67
68 // Tecplot header info
69 outfile << "ZONE I=" << n_plot << std::endl;
70
71 // Find the dimension of the nodes
72 unsigned n_dim = this->nodal_dimension();
73
74 // Loop over plot points
75 for (unsigned l = 0; l < n_plot; l++)
76 {
77 s[0] = l * 1.0 / (n_plot - 1);
78 // Output the coordinates
79 for (unsigned i = 0; i < n_dim; i++)
80 {
81 outfile << interpolated_x(s, i) << " ";
82 }
83 outfile << std::endl;
84 }
85 outfile << std::endl;
86 }
87
88
89 //=======================================================================
90 /// C style output function for general 1D TElements
91 //=======================================================================
92 template<unsigned NNODE_1D>
94 {
95 output(file_pt, NNODE_1D);
96 }
97
98 //=======================================================================
99 /// C style output function for n_plot points in each coordinate direction
100 //=======================================================================
101 template<unsigned NNODE_1D>
102 void TElement<1, NNODE_1D>::output(FILE* file_pt, const unsigned& n_plot)
103 {
104 // Local variables
105 Vector<double> s(1);
106
107 // Tecplot header info
108 // outfile << "ZONE I=" << n_plot << std::endl;
109 fprintf(file_pt, "ZONE I=%i\n", n_plot);
110
111 // Find the dimension of the first node
112 unsigned n_dim = this->nodal_dimension();
113
114 // Loop over plot points
115 for (unsigned l = 0; l < n_plot; l++)
116 {
117 s[0] = l * 1.0 / (n_plot - 1);
118
119 // Output the coordinates
120 for (unsigned i = 0; i < n_dim; i++)
121 {
122 // outfile << interpolated_x(s,i) << " " ;
123 fprintf(file_pt, "%g ", interpolated_x(s, i));
124 }
125 // outfile << std::endl;
126 fprintf(file_pt, "\n");
127 }
128 // outfile << std::endl;
129 fprintf(file_pt, "\n");
130 }
131
132 //=============================================================
133 /// Namespace for helper functions that return the local
134 /// coordinates in the bulk elements
135 //==============================================================
136 namespace TElement1FaceToBulkCoordinates
137 {
138 /// The translation scheme for the face s0 = 0.0
139 void face0(const Vector<double>& s, Vector<double>& s_bulk)
140 {
141 s_bulk[0] = 0.0;
142 }
143
144 /// The translation scheme for the face s0 = 1.0
145 void face1(const Vector<double>& s, Vector<double>& s_bulk)
146 {
147 s_bulk[0] = 1.0;
148 }
149 } // namespace TElement1FaceToBulkCoordinates
150
151
152 //=============================================================
153 /// Namespace for helper functions that calculate derivatives
154 /// of the local coordinates in the bulk elements wrt the
155 /// local coordinates in the face element.
156 //=============================================================
157 namespace TElement1BulkCoordinateDerivatives
158 {
159 /// Function for both faces -- the bulk coordinate is fixed on both
161 DenseMatrix<double>& dsbulk_dsface,
162 unsigned& interior_direction)
163 {
164 // Bulk coordinate s[0] does not vary along the face
165 dsbulk_dsface(0, 0) = 0.0;
166
167 // The interior direction is given by s[0]
168 interior_direction = 0;
169 }
170 } // namespace TElement1BulkCoordinateDerivatives
171
172 //===========================================================
173 /// Function to setup geometrical information for lower-dimensional
174 /// FaceElements (which are Point Elements).
175 //===========================================================
176 template<unsigned NNODE_1D>
178 FaceElement* face_element_pt)
179 {
180 // Overload the nodal dimension
181 face_element_pt->set_nodal_dimension(nodal_dimension());
182
183 // Set the pointer to the "bulk" element
184 face_element_pt->bulk_element_pt() = this;
185
186#ifdef OOMPH_HAS_MPI
187 // Pass on non-halo proc ID
188 face_element_pt->set_halo(Non_halo_proc_ID);
189#endif
190
191 // Resize the storage for the original number of values at the (one and
192 // only) node of the face element.
193 face_element_pt->nbulk_value_resize(1);
194
195 // Resize the storage for the bulk node number corresponding to the (one
196 // and only) node of the face element
197 face_element_pt->bulk_node_number_resize(1);
198
199 // Set the face index in the face element
200 face_element_pt->face_index() = face_index;
201
202 // Now set up the node pointer
203 // The convention is that the "normal", should always point
204 // out of the element
205 switch (face_index)
206 {
207 // Bottom, normal sign is negative (coordinate points into element)
208 case (-1):
209 face_element_pt->node_pt(0) = node_pt(0);
210 face_element_pt->bulk_node_number(0) = 0;
211 face_element_pt->normal_sign() = -1;
212
213 // Set the pointer to the function that determines the bulk coordinates
214 // in the face element
215 face_element_pt->face_to_bulk_coordinate_fct_pt() =
217
218 // Set the pointer to the function that determines the mapping of
219 // derivatives
220 face_element_pt->bulk_coordinate_derivatives_fct_pt() =
222
223 // Set the number of values stored when the node is part of the "bulk"
224 // element. The required_nvalue() must be used, rather than nvalue(),
225 // because otherwise nodes on boundaries will be resized multiple
226 // times. If you want any other behaviour, you MUST set nbulk_value()
227 // manually after construction of your specific face element.
228 face_element_pt->nbulk_value(0) = required_nvalue(0);
229 break;
230
231 // Top, normal sign is positive (coordinate points out of element)
232 case (1):
233 face_element_pt->node_pt(0) = node_pt(NNODE_1D - 1);
234 face_element_pt->bulk_node_number(0) = NNODE_1D - 1;
235 face_element_pt->normal_sign() = +1;
236
237 // Set the pointer to the function that determines the bulk coordinates
238 // in the face element
239 face_element_pt->face_to_bulk_coordinate_fct_pt() =
241
242 // Set the pointer to the function that determines the mapping of
243 // derivatives
244 face_element_pt->bulk_coordinate_derivatives_fct_pt() =
246
247
248 // Set the number of values stored when the node is part of the "bulk"
249 // element.
250 face_element_pt->nbulk_value(0) = required_nvalue(NNODE_1D - 1);
251 break;
252
253 // All other cases throw an error
254 default:
255 std::ostringstream error_message;
256 error_message << "Face_index should only take "
257 << "the values +/-1, not " << face_index << std::endl;
258
259 throw OomphLibError(error_message.str(),
260 OOMPH_CURRENT_FUNCTION,
261 OOMPH_EXCEPTION_LOCATION);
262 }
263 }
264
265
266 /// /////////////////////////////////////////////////////////////
267 // 2D Telements
268 /// /////////////////////////////////////////////////////////////
269
270 /// Assign the nodal translation schemes
271 template<>
272 const unsigned TElement<2, 2>::Node_on_face[3][2] = {{2, 1}, {2, 0}, {0, 1}};
273
274 template<>
275 const unsigned TElement<2, 3>::Node_on_face[3][3] = {
276 {2, 4, 1}, {2, 5, 0}, {0, 3, 1}};
277
278 template<>
279 const unsigned TElement<2, 4>::Node_on_face[3][4] = {
280 {2, 6, 5, 1}, {2, 7, 8, 0}, {0, 3, 4, 1}};
281
282
283 //===================================================================
284 /// Namespace for the functions that translate local face coordinates
285 /// to the coordinates in the bulk element
286 //==================================================================
287 namespace TElement2FaceToBulkCoordinates
288 {
289 /// The translation scheme for the face s0 = 0
290 void face0(const Vector<double>& s, Vector<double>& s_bulk)
291 {
292 s_bulk[0] = 0.0;
293 s_bulk[1] = s[0];
294 }
295
296 /// The translation scheme for the face s1 = 0
297 void face1(const Vector<double>& s, Vector<double>& s_bulk)
298 {
299 s_bulk[0] = s[0];
300 s_bulk[1] = 0.0;
301 }
302
303 /// The translation scheme for the face s2 = 0
304 void face2(const Vector<double>& s, Vector<double>& s_bulk)
305 {
306 s_bulk[0] = 1.0 - s[0];
307 s_bulk[1] = s[0];
308 }
309 } // namespace TElement2FaceToBulkCoordinates
310
311
312 //=============================================================
313 /// Namespace for helper functions that calculate derivatives
314 /// of the local coordinates in the bulk elements wrt the
315 /// local coordinates in the face element.
316 //=============================================================
317 namespace TElement2BulkCoordinateDerivatives
318 {
319 /// Function for the "left" face along which s0 is fixed
320 void face0(const Vector<double>& s,
321 DenseMatrix<double>& dsbulk_dsface,
322 unsigned& interior_direction)
323 {
324 // Bulk coordinate s[0] does not vary along the face
325 dsbulk_dsface(0, 0) = 0.0;
326 // Bulk coordinate s[1] is the face coordinate
327 dsbulk_dsface(1, 0) = 1.0;
328
329 // The interior direction is given by s[0]
330 interior_direction = 0;
331 }
332
333
334 /// Function for the "bottom" face along which s1 is fixed
335 void face1(const Vector<double>& s,
336 DenseMatrix<double>& dsbulk_dsface,
337 unsigned& interior_direction)
338 {
339 // Bulk coordinate s[0] is the face coordinate
340 dsbulk_dsface(0, 0) = 1.0;
341 // Bulk coordinate s[1] does not vary along the face
342 dsbulk_dsface(1, 0) = 0.0;
343
344 // The interior direction is given by s[1]
345 interior_direction = 1;
346 }
347
348 /// Function for the sloping face
349 void face2(const Vector<double>& s,
350 DenseMatrix<double>& dsbulk_dsface,
351 unsigned& interior_direction)
352 {
353 // Bulk coordinate s[0] decreases along the face
354 dsbulk_dsface(0, 0) = -1.0;
355 // Bulk coordinate s[1] increases along the face
356 dsbulk_dsface(1, 0) = 1.0;
357
358 // The interior direction is given by s[0] (or s[1])
359 interior_direction = 0;
360 }
361
362 } // namespace TElement2BulkCoordinateDerivatives
363
364
365 //=======================================================================
366 /// Function to setup geometrical information for lower-dimensional
367 /// FaceElements (which are of type TElement<2,NNODE_1D>).
368 //=======================================================================
369 template<unsigned NNODE_1D>
371 FaceElement* face_element_pt)
372 {
373 // Set the nodal dimension from the first node
374 face_element_pt->set_nodal_dimension(nodal_dimension());
375
376 // Set the pointer to the orginal "bulk" element
377 face_element_pt->bulk_element_pt() = this;
378
379#ifdef OOMPH_HAS_MPI
380 // Pass on non-halo proc ID
381 face_element_pt->set_halo(Non_halo_proc_ID);
382#endif
383
384 // Calculate the number of nodes in the face element
385 const unsigned n_face_nodes = NNODE_1D;
386
387 // Resize storage for the number of values originally stored
388 // at the face element's nodes.
389 face_element_pt->nbulk_value_resize(n_face_nodes);
390
391 // Resize storage for the bulk node numbers corresponding to
392 // the nodes of the face
393 face_element_pt->bulk_node_number_resize(n_face_nodes);
394
395 // Set the face index in the face element
396 face_element_pt->face_index() = face_index;
397
398 // So the faces are
399 // 0 : s_0 fixed
400 // 1 : s_1 fixed
401 // 2 : sloping face
402
403 // Copy nodes across to the face
404 for (unsigned i = 0; i < n_face_nodes; i++)
405 {
406 // The number is just offset by one
407 unsigned bulk_number = Node_on_face[face_index][i];
408 face_element_pt->node_pt(i) = node_pt(bulk_number);
409 face_element_pt->bulk_node_number(i) = bulk_number;
410 // set the number of values originally stored at this node
411 face_element_pt->nbulk_value(i) = required_nvalue(bulk_number);
412 }
413
414 // Now set up the node pointers and the normal vectors
415 switch (face_index)
416 {
417 //
418 //-----The face s0 is constant
419 case 0:
420
421 // Set the pointer to the function that determines the bulk
422 // coordinates in the face element
423 face_element_pt->face_to_bulk_coordinate_fct_pt() =
425
426 // Set the pointer to the function that determines the mapping of
427 // derivatives
428 face_element_pt->bulk_coordinate_derivatives_fct_pt() =
430
431 // Outer unit normal is the positive cross product of two in plane
432 // tangent vectors
433 face_element_pt->normal_sign() = +1;
434
435 break;
436
437 //-----The face s1 is constant
438 case 1:
439
440 // Set the pointer to the function that determines the bulk
441 // coordinates in the face element
442 face_element_pt->face_to_bulk_coordinate_fct_pt() =
444
445 // Set the pointer to the function that determines the mapping of
446 // derivatives
447 face_element_pt->bulk_coordinate_derivatives_fct_pt() =
449
450 // Outer unit normal is the positive cross product of two in plane
451 // tangent vectors
452 face_element_pt->normal_sign() = +1;
453
454 break;
455
456 //
457 //-----The face s2 is constant
458 case 2:
459
460 // Set the pointer to the function that determines the bulk
461 // coordinates in the face element
462 face_element_pt->face_to_bulk_coordinate_fct_pt() =
464
465 // Set the pointer to the function that determines the mapping of
466 // derivatives
467 face_element_pt->bulk_coordinate_derivatives_fct_pt() =
469
470 // Outer unit normal is the negative cross product of two in plane
471 // tangent vectors
472 face_element_pt->normal_sign() = -1;
473
474 break;
475
476 default:
477
478 std::ostringstream error_message;
479 error_message << "Face_index should only take "
480 << "the values 0, 1 or 2 not " << face_index << std::endl;
481
482 throw OomphLibError(error_message.str(),
483 OOMPH_CURRENT_FUNCTION,
484 OOMPH_EXCEPTION_LOCATION);
485 } // end switch
486 }
487
488
489 //=======================================================================
490 /// The output function for TElement<2,NNODE_1D>
491 //=======================================================================
492 template<unsigned NNODE_1D>
493 void TElement<2, NNODE_1D>::output(std::ostream& outfile)
494 {
495 // QUEHACERES want to perform same output, but at each node
496 output(outfile, NNODE_1D);
497 }
498
499
500 //=======================================================================
501 /// The output function for TElement<2,NNODE_1D>
502 //=======================================================================
503 template<unsigned NNODE_1D>
504 void TElement<2, NNODE_1D>::output(std::ostream& outfile,
505 const unsigned& nplot)
506 {
507 // Vector of local coordinates
508 Vector<double> s(2);
509
510 // Get the dimension of the node
511 unsigned n_dim = this->nodal_dimension();
512
513 // Tecplot header info
514 outfile << tecplot_zone_string(nplot);
515
516 // Loop over plot points
517 unsigned num_plot_points = nplot_points(nplot);
518 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
519 {
520 // Get local coordinates of plot point
521 get_s_plot(iplot, nplot, s);
522
523 for (unsigned i = 0; i < n_dim; i++)
524 {
525 outfile << interpolated_x(s, i) << " ";
526 }
527 outfile << std::endl;
528 }
529
530 // Write tecplot footer (e.g. FE connectivity lists)
531 write_tecplot_zone_footer(outfile, nplot);
532 }
533
534 //=======================================================================
535 /// The C-style output function for TElement<2,NNODE_1D>
536 //=======================================================================
537 template<unsigned NNODE_1D>
539 {
540 output(file_pt, NNODE_1D);
541 }
542
543
544 //=======================================================================
545 /// The C-style output function for TElement<2,NNODE_1D>
546 //=======================================================================
547 template<unsigned NNODE_1D>
548 void TElement<2, NNODE_1D>::output(FILE* file_pt, const unsigned& nplot)
549 {
550 // Vector of local coordinates
551 Vector<double> s(2);
552
553 // Find the dimensions of the nodes
554 unsigned n_dim = this->nodal_dimension();
555
556 // Tecplot header info
557 fprintf(file_pt, "%s \n", tecplot_zone_string(nplot).c_str());
558
559 // Loop over plot points
560 unsigned num_plot_points = nplot_points(nplot);
561 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
562 {
563 // Get local coordinates of plot point
564 get_s_plot(iplot, nplot, s);
565
566 for (unsigned i = 0; i < n_dim; i++)
567 {
568 fprintf(file_pt, "%g ", interpolated_x(s, i));
569 // outfile << interpolated_x(s,i) << " ";
570 }
571 fprintf(file_pt, "\n");
572 // outfile << std::endl;
573 }
574
575 // Write tecplot footer (e.g. FE connectivity lists)
576 write_tecplot_zone_footer(file_pt, nplot);
577 }
578
579
580 /// ////////////////////////////////////////////////////////////////////
581 /// ////////////////////////////////////////////////////////////////////
582 /// ////////////////////////////////////////////////////////////////////
583
584
585 //=======================================================================
586 /// The output function for TElement<3,NNODE_1D>
587 //=======================================================================
588 template<unsigned NNODE_1D>
589 void TElement<3, NNODE_1D>::output(std::ostream& outfile)
590 {
591 output(outfile, NNODE_1D);
592 }
593
594
595 //=======================================================================
596 /// The output function for TElement<3,NNODE_1D>
597 //=======================================================================
598 template<unsigned NNODE_1D>
599 void TElement<3, NNODE_1D>::output(std::ostream& outfile,
600 const unsigned& nplot)
601 {
602 // Vector of local coordinates
603 Vector<double> s(3);
604
605 // Find the dimension of the nodes
606 unsigned n_dim = this->nodal_dimension();
607
608 // Tecplot header info
609 outfile << tecplot_zone_string(nplot);
610
611 // Loop over plot points
612 unsigned num_plot_points = nplot_points(nplot);
613 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
614 {
615 // Get local coordinates of plot point
616 get_s_plot(iplot, nplot, s);
617
618 for (unsigned i = 0; i < n_dim; i++)
619 {
620 outfile << interpolated_x(s, i) << " ";
621 }
622 outfile << "\n";
623 }
624
625 // Write tecplot footer (e.g. FE connectivity lists)
626 write_tecplot_zone_footer(outfile, nplot);
627 }
628
629
630 //=======================================================================
631 /// The C-style output function for TElement<3,NNODE_1D>
632 //=======================================================================
633 template<unsigned NNODE_1D>
635 {
636 output(file_pt, NNODE_1D);
637 }
638
639
640 //=======================================================================
641 /// The C-style output function for TElement<3,NNODE_1D>
642 //=======================================================================
643 template<unsigned NNODE_1D>
644 void TElement<3, NNODE_1D>::output(FILE* file_pt, const unsigned& nplot)
645 {
646 // Vector of local coordinates
647 Vector<double> s(3);
648
649 // Find the dimension of the nodes
650 unsigned n_dim = this->nodal_dimension();
651
652 // Tecplot header info
653 fprintf(file_pt, "%s \n", tecplot_zone_string(nplot).c_str());
654
655 // Loop over plot points
656 unsigned num_plot_points = nplot_points(nplot);
657 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
658 {
659 // Get local coordinates of plot point
660 get_s_plot(iplot, nplot, s);
661
662 for (unsigned i = 0; i < n_dim; i++)
663 {
664 fprintf(file_pt, "%g ", interpolated_x(s, i));
665 // outfile << interpolated_x(s,i) << " ";
666 }
667 fprintf(file_pt, "\n");
668 // outfile << std::endl;
669 }
670
671 // Write tecplot footer (e.g. FE connectivity lists)
672 write_tecplot_zone_footer(file_pt, nplot);
673 }
674
675 //===================================================================
676 /// Namespace for the functions that translate local face coordinates
677 /// to the coordinates in the bulk element
678 //==================================================================
679 namespace TElement3FaceToBulkCoordinates
680 {
681 /// The translation scheme for the face s0 = 0
682 void face0(const Vector<double>& s, Vector<double>& s_bulk)
683 {
684 s_bulk[0] = 0.0;
685 s_bulk[1] = s[0];
686 s_bulk[2] = s[1];
687 }
688
689 /// The translation scheme for the face s1 = 0
690 void face1(const Vector<double>& s, Vector<double>& s_bulk)
691 {
692 s_bulk[0] = s[0];
693 s_bulk[1] = 0.0;
694 s_bulk[2] = s[1];
695 }
696
697 /// The translation scheme for the face s2 = 0
698 void face2(const Vector<double>& s, Vector<double>& s_bulk)
699 {
700 s_bulk[0] = s[0];
701 s_bulk[1] = s[1];
702 s_bulk[2] = 0.0;
703 }
704
705 /// The translation scheme for the sloping face
706 void face3(const Vector<double>& s, Vector<double>& s_bulk)
707 {
708 s_bulk[0] = 1 - s[0] - s[1];
709 s_bulk[1] = s[0];
710 s_bulk[2] = s[1];
711 }
712
713 } // namespace TElement3FaceToBulkCoordinates
714
715 /// Assign the nodal translation scheme for linear elements
716 template<>
717 const unsigned TElement<3, 2>::Node_on_face[4][3] = {
718 {1, 2, 3}, {0, 2, 3}, {0, 1, 3}, {1, 2, 0}};
719
720 /// Assign the nodal translation scheme for quadratic elements
721 template<>
722 const unsigned TElement<3, 3>::Node_on_face[4][6] = {{1, 2, 3, 7, 8, 9},
723 {0, 2, 3, 5, 8, 6},
724 {0, 1, 3, 4, 9, 6},
725 {1, 2, 0, 7, 5, 4}};
726
727
728 //=======================================================================
729 /// Function to setup geometrical information for lower-dimensional
730 /// FaceElements (which are of type TElement<2,NNODE_1D>).
731 //=======================================================================
732 template<unsigned NNODE_1D>
734 FaceElement* face_element_pt)
735 {
736 // Set the nodal dimension
737 face_element_pt->set_nodal_dimension(nodal_dimension());
738
739 // Set the pointer to the orginal "bulk" element
740 face_element_pt->bulk_element_pt() = this;
741
742#ifdef OOMPH_HAS_MPI
743 // Pass on non-halo proc ID
744 face_element_pt->set_halo(Non_halo_proc_ID);
745#endif
746
747 // Calculate the number of nodes in the face element
748 const unsigned n_face_nodes = (NNODE_1D * (NNODE_1D + 1)) / 2;
749
750 // Resize storage for the number of values originally stored
751 // at the face element's nodes.
752 face_element_pt->nbulk_value_resize(n_face_nodes);
753
754 // Resize storage for the bulk node numbers corresponding to
755 // the nodes of the face
756 face_element_pt->bulk_node_number_resize(n_face_nodes);
757
758 // Set the face index in the face element
759 face_element_pt->face_index() = face_index;
760
761 // So the faces are
762 // 0 : s_0 fixed
763 // 1 : s_1 fixed
764 // 2 : s_2 fixed
765 // 3 : sloping face
766
767 // Copy nodes across to the face
768 for (unsigned i = 0; i < n_face_nodes; i++)
769 {
770 // The number is just offset by one
771 unsigned bulk_number = Node_on_face[face_index][i];
772 face_element_pt->node_pt(i) = node_pt(bulk_number);
773 face_element_pt->bulk_node_number(i) = bulk_number;
774 // set the number of values originally stored at this node
775 face_element_pt->nbulk_value(i) = required_nvalue(bulk_number);
776 }
777
778 // Now set up the node pointers and the normal vectors
779 switch (face_index)
780 {
781 //
782 //-----The face s0 is constant
783 case 0:
784
785 // Set the pointer to the function that determines the bulk
786 // coordinates in the face element
787 face_element_pt->face_to_bulk_coordinate_fct_pt() =
789
790 // Outer unit normal is the negative cross product of two in plane
791 // tangent vectors
792 face_element_pt->normal_sign() = -1;
793
794 break;
795
796 //-----The face s1 is constant
797 case 1:
798
799 // Set the pointer to the function that determines the bulk
800 // coordinates in the face element
801 face_element_pt->face_to_bulk_coordinate_fct_pt() =
803
804 // Outer unit normal is the positive cross product of two in plane
805 // tangent vectors
806 face_element_pt->normal_sign() = +1;
807
808 break;
809
810 //
811 //-----The face s2 is constant
812 case 2:
813
814 // Set the pointer to the function that determines the bulk
815 // coordinates in the face element
816 face_element_pt->face_to_bulk_coordinate_fct_pt() =
818
819 // Outer unit normal is the negative cross product of two in plane
820 // tangent vectors
821 face_element_pt->normal_sign() = -1;
822
823 break;
824
825 //-----The sloping face of the tetrahedron
826 case 3:
827
828 // Set the pointer to the function that determines the bulk
829 // coordinates in the face element
830 face_element_pt->face_to_bulk_coordinate_fct_pt() =
832
833 // Outer unit normal is the positive cross product of two in plane
834 // tangent vectors
835 face_element_pt->normal_sign() = +1;
836
837 break;
838
839
840 default:
841
842 std::ostringstream error_message;
843 error_message << "Face_index should only take "
844 << "the values 0, 1, 2 or 3, not " << face_index
845 << std::endl;
846
847 throw OomphLibError(error_message.str(),
848 OOMPH_CURRENT_FUNCTION,
849 OOMPH_EXCEPTION_LOCATION);
850 } // end switch
851 }
852
853
854 //==================================================================
855 // Default integration scheme for the TBubbleEnrichedElement<2,3>
856 //===================================================================
857 template<unsigned DIM>
860
861 //===================================================================
862 // Central node on the face of the TBubbleEnrichedelement<2,3>
863 //===================================================================
864 template<>
866 4, 5, 3};
867
868
869 //================================================================
870 /// The face element for is the same in the two-dimesional case
871 //================================================================
872 template<>
874 const int& face_index, FaceElement* face_element_pt)
875 {
876 TElement<2, 3>::build_face_element(face_index, face_element_pt);
877 }
878
879
880 //===================================================================
881 // Central node on the face of the TBubbleEnrichedelement<3,3>
882 //===================================================================
883 template<>
885 13, 12, 10, 11};
886
887
888 //=======================================================================
889 /// Function to setup geometrical information for lower-dimensional
890 /// FaceElements (which are of type TBubbleEnrichedElement<2,3>).
891 //=======================================================================
892 template<>
894 const int& face_index, FaceElement* face_element_pt)
895 {
896 // Call the standard unenriched build function
897 TElement<3, 3>::build_face_element(face_index, face_element_pt);
898
899 // Set the enriched number of total face nodes
900 const unsigned n_face_nodes = 7;
901
902 // Resize storage for the number of values originally stored
903 // at the face element's nodes.
904 face_element_pt->nbulk_value_resize(n_face_nodes);
905
906 // Resize storage for the bulk node numbers corresponding to
907 // the nodes of the face
908 face_element_pt->bulk_node_number_resize(n_face_nodes);
909
910 // So the faces are
911 // 0 : s_0 fixed
912 // 1 : s_1 fixed
913 // 2 : s_2 fixed
914 // 3 : sloping face
915
916 // Copy central node across
917 unsigned bulk_number = Central_node_on_face[face_index];
918 face_element_pt->node_pt(n_face_nodes - 1) = node_pt(bulk_number);
919 face_element_pt->bulk_node_number(n_face_nodes - 1) = bulk_number;
920 // set the number of values originally stored at this node
921 face_element_pt->nbulk_value(n_face_nodes - 1) =
922 required_nvalue(bulk_number);
923 }
924
925 /// //////////////////////////////////////////////////////////////////////
926 /// //////////////////////////////////////////////////////////////////////
927 /// //////////////////////////////////////////////////////////////////////
928
929 //===========================================================
930 /// Final override for
931 /// function to setup geometrical information for lower-dimensional
932 /// FaceElements (which are of type SolidTBubbleEnrichedElement<1,3>).
933 //===========================================================
934 template<>
936 const int& face_index, FaceElement* face_element_pt)
937 {
938 // Build the standard non-solid FaceElement
940 face_element_pt);
941
942 // Set the Lagrangian dimension from the first node of the present element
943 dynamic_cast<SolidFiniteElement*>(face_element_pt)
944 ->set_lagrangian_dimension(
945 static_cast<SolidNode*>(node_pt(0))->nlagrangian());
946 }
947
948
949 //===========================================================
950 /// Final override for
951 /// function to setup geometrical information for lower-dimensional
952 /// FaceElements (which are of type SolidTBubbleEnrichedElement<2,3>).
953 //===========================================================
954 template<>
956 const int& face_index, FaceElement* face_element_pt)
957 {
958 // Build the standard non-solid FaceElement
960 face_element_pt);
961
962 // Set the Lagrangian dimension from the first node of the present element
963 dynamic_cast<SolidFiniteElement*>(face_element_pt)
964 ->set_lagrangian_dimension(
965 static_cast<SolidNode*>(node_pt(0))->nlagrangian());
966 }
967
968
969 //===================================================================
970 // Build required templates
971 //===================================================================
972 template class TElement<1, 2>;
973 template class TElement<1, 3>;
974 template class TElement<1, 4>;
975 template class TElement<2, 2>;
976 template class TElement<2, 3>;
977 template class TElement<2, 4>;
978 template class TElement<3, 2>;
979 template class TElement<3, 3>;
980 template class TBubbleEnrichedElement<2, 3>;
981 template class TBubbleEnrichedElement<3, 3>;
984 template class TBubbleEnrichedGauss<2, 3>;
985 template class TBubbleEnrichedGauss<3, 3>;
986} // namespace oomph
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
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
FiniteElement *& bulk_element_pt()
Pointer to higher-dimensional "bulk" element.
Definition: elements.h:4735
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....
////////////////////////////////////////////////////////////////////// //////////////////////////////...
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
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: Telements.h:3926
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: Telements.h:3570
void build_face_element(const int &face_index, FaceElement *face_element_pt)
The face element for is the same in the two-dimesional case.
Definition: Telements.cc:873
Specialisation for two-dimensional elements, in which the highest order polynomial is cubic,...
Definition: Telements.h:3591
Specialisation for three-dimensional elements, in which the highest order polynomial is quartic,...
Definition: Telements.h:3603
Define integration schemes that are required to exactly integrate the mass matrices of the bubble-enr...
Definition: Telements.h:3581
static TGauss< 1, NNODE_1D > Default_integration_scheme
Default integration rule: Gaussian integration of same 'order' as the element.
Definition: Telements.h:1228
static TGauss< 2, NNODE_1D > Default_integration_scheme
Default integration rule: Gaussian integration of same 'order' as the element.
Definition: Telements.h:1515
static TGauss< 3, NNODE_1D > Default_integration_scheme
Default integration rule: Gaussian integration of same 'order' as the element.
Definition: Telements.h:3052
General TElement class.
Definition: Telements.h:1208
const unsigned Node_on_face[3][2]
Assign the nodal translation schemes.
Definition: Telements.cc:272
void output()
Doc the command line arguments.
void faces0(const Vector< double > &s, DenseMatrix< double > &dsbulk_dsface, unsigned &interior_direction)
Function for both faces – the bulk coordinate is fixed on both.
Definition: Telements.cc:160
void face0(const Vector< double > &s, Vector< double > &s_bulk)
The translation scheme for the face s0 = 0.0.
Definition: Telements.cc:139
void face1(const Vector< double > &s, Vector< double > &s_bulk)
The translation scheme for the face s0 = 1.0.
Definition: Telements.cc:145
void face0(const Vector< double > &s, DenseMatrix< double > &dsbulk_dsface, unsigned &interior_direction)
Function for the "left" face along which s0 is fixed.
Definition: Telements.cc:320
void face1(const Vector< double > &s, DenseMatrix< double > &dsbulk_dsface, unsigned &interior_direction)
Function for the "bottom" face along which s1 is fixed.
Definition: Telements.cc:335
void face2(const Vector< double > &s, DenseMatrix< double > &dsbulk_dsface, unsigned &interior_direction)
Function for the sloping face.
Definition: Telements.cc:349
void face1(const Vector< double > &s, Vector< double > &s_bulk)
The translation scheme for the face s1 = 0.
Definition: Telements.cc:297
void face2(const Vector< double > &s, Vector< double > &s_bulk)
The translation scheme for the face s2 = 0.
Definition: Telements.cc:304
void face0(const Vector< double > &s, Vector< double > &s_bulk)
The translation scheme for the face s0 = 0.
Definition: Telements.cc:290
void face1(const Vector< double > &s, Vector< double > &s_bulk)
The translation scheme for the face s1 = 0.
Definition: Telements.cc:690
void face2(const Vector< double > &s, Vector< double > &s_bulk)
The translation scheme for the face s2 = 0.
Definition: Telements.cc:698
void face0(const Vector< double > &s, Vector< double > &s_bulk)
The translation scheme for the face s0 = 0.
Definition: Telements.cc:682
void face3(const Vector< double > &s, Vector< double > &s_bulk)
The translation scheme for the sloping face.
Definition: Telements.cc:706
//////////////////////////////////////////////////////////////////// ////////////////////////////////...