Qelements.h
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// Header functions for classes that define Qelements
27// Include guards to prevent multiple inclusion of the header
28#ifndef OOMPH_QELEMENT_HEADER
29#define OOMPH_QELEMENT_HEADER
30
31// Config header generated by autoconfig
32#ifdef HAVE_CONFIG_H
33#include <oomph-lib-config.h>
34#endif
35
36
37#ifdef OOMPH_HAS_MPI
38#include "mpi.h"
39#endif
40
41// oomph-lib headers
42#include "Vector.h"
43#include "shape.h"
44#include "integral.h"
45#include "timesteppers.h"
46#include "elements.h"
47#include "macro_element.h"
48
50
51
52namespace oomph
53{
54 /// ///////////////////////////////////////////////////////////////////
55 /// ///////////////////////////////////////////////////////////////////
56 /// ///////////////////////////////////////////////////////////////////
57
58 //========================================================================
59 /// Empty base class for Qelements (created so that
60 /// we can use dynamic_cast<>() to figure out if a an element
61 /// is a Qelement (from a purely geometric point of view).
62 //========================================================================
63 class QElementGeometricBase : public virtual FiniteElement
64 {
65 public:
66 /// Empty default constructor
68
69 /// Broken copy constructor
71
72 /// Broken assignment operator
73 // Commented out broken assignment operator because this can lead to a
74 // conflict warning when used in the virtual inheritence hierarchy.
75 // Essentially the compiler doesn't realise that two separate
76 // implementations of the broken function are the same and so, quite
77 // rightly, it shouts.
78 /*void operator=(const QElementGeometricBase&) = delete;*/
79 };
80
81
82 /// ///////////////////////////////////////////////////////////////////
83 /// ///////////////////////////////////////////////////////////////////
84 /// ///////////////////////////////////////////////////////////////////
85
86
87 //========================================================================
88 /// Base class for Qelements
89 //========================================================================
90 class QElementBase : public virtual QElementGeometricBase
91 {
92 public:
93 /// Constructor: Initialise pointers to macro element reference coords
95
96 /// Broken copy constructor
97 QElementBase(const QElementBase&) = delete;
98
99 /// Broken assignment operator
100 /*void operator=(const QElementBase&) = delete;*/
101
102 /// Destructor: Kill storage for macro element reference coords
104 {
105 // Can be deleted blindly as they were nulled initially
106 delete S_macro_ll_pt;
107 S_macro_ll_pt = 0;
108 delete S_macro_ur_pt;
109 S_macro_ur_pt = 0;
110 }
111
112 /// Check whether the local coordinate are valid or not
114 {
115 unsigned ncoord = dim();
116 for (unsigned i = 0; i < ncoord; i++)
117 {
118 // We're outside
119 if ((s[i] - s_max() > 0.0) || (s_min() - s[i] > 0.0))
120 {
121 return false;
122 }
123 }
124 return true;
125 }
126
127 /// Adjust local coordinates so that they're located inside
128 /// the element
130 {
131 unsigned ncoord = dim();
132 for (unsigned i = 0; i < ncoord; i++)
133 {
134 // Adjust to move it onto the boundary
135 if (s[i] > s_max()) s[i] = s_max();
136 if (s[i] < s_min()) s[i] = s_min();
137 }
138 }
139
140
141 /// Set pointer to macro element also sets up storage for the
142 /// reference coordinates and initialises them
144 {
145 // Get the spatial dimension (= number of local and macro element coords)
146 unsigned n_dim = dim();
147
148 // Create storage if none has been allocated
149 if (S_macro_ll_pt == 0)
150 {
151 S_macro_ll_pt = new Vector<double>(n_dim);
152 }
153 // Otherwise resize the allocated storage
154 else
155 {
156 S_macro_ll_pt->resize(n_dim);
157 }
158
159 // Create storage if none has been allocated
160 if (S_macro_ur_pt == 0)
161 {
162 S_macro_ur_pt = new Vector<double>(n_dim);
163 }
164 // Otherwise resize the allocated storage
165 else
166 {
167 S_macro_ur_pt->resize(n_dim);
168 }
169
170 // Initialise the vertex coordinates in the macro element
171 // Default: The element is unrefined and hence its vertices are those
172 // of the macro element itself
173 for (unsigned i = 0; i < n_dim; i++)
174 {
175 s_macro_ll(i) = -1.0;
176 s_macro_ur(i) = 1.0;
177 }
178
179 /// Call the corresponding function in the FiniteElement base class
181 }
182
183
184 /// Access fct to the i-th coordinate of the element's
185 /// "lower left" vertex in the associated MacroElement
186 double& s_macro_ll(const unsigned& i)
187 {
188#ifdef PARANOID
189 if (S_macro_ll_pt == 0)
190 {
191 throw OomphLibError("S_macro_ll_pt has not been set\n",
192 OOMPH_CURRENT_FUNCTION,
193 OOMPH_EXCEPTION_LOCATION);
194 }
195#endif
196 return (*S_macro_ll_pt)[i];
197 }
198
199
200 /// Access fct to the i-th coordinate of the element's
201 /// "upper right" vertex in the associated MacroElement
202 double& s_macro_ur(const unsigned& i)
203 {
204#ifdef PARANOID
205 if (S_macro_ur_pt == 0)
206 {
207 throw OomphLibError("S_macro_ur_pt has not been set\n",
208 OOMPH_CURRENT_FUNCTION,
209 OOMPH_EXCEPTION_LOCATION);
210 }
211#endif
212 return (*S_macro_ur_pt)[i];
213 }
214
215
216 /// Access fct to the i-th coordinate of the element's
217 /// "lower left" vertex in the associated MacroElement. (const version)
218 double s_macro_ll(const unsigned& i) const
219 {
220#ifdef PARANOID
221 if (S_macro_ll_pt == 0)
222 {
223 throw OomphLibError("S_macro_ll_pt has not been set\n",
224 OOMPH_CURRENT_FUNCTION,
225 OOMPH_EXCEPTION_LOCATION);
226 }
227#endif
228 return (*S_macro_ll_pt)[i];
229 }
230
231
232 /// Access fct to the i-th coordinate of the element's
233 /// "upper right" vertex in the associated MacroElement. (const version)
234 double s_macro_ur(const unsigned& i) const
235 {
236#ifdef PARANOID
237 if (S_macro_ur_pt == 0)
238 {
239 throw OomphLibError("S_macro_ur_pt has not been set\n",
240 OOMPH_CURRENT_FUNCTION,
241 OOMPH_EXCEPTION_LOCATION);
242 }
243#endif
244 return (*S_macro_ur_pt)[i];
245 }
246
247 /// Global coordinates as function of local coordinates.
248 /// using the macro element representation
250 Vector<double>& x) const
251 {
252 // Check that there is a macro element
253 if (Macro_elem_pt == 0)
254 {
255 throw OomphLibError("Macro Element pointer not set in this element\n",
256 OOMPH_CURRENT_FUNCTION,
257 OOMPH_EXCEPTION_LOCATION);
258 }
259
260 // Use macro element representation
261 unsigned el_dim = dim();
262 Vector<double> s_macro(el_dim);
263 for (unsigned i = 0; i < el_dim; i++)
264 {
265 s_macro[i] =
266 s_macro_ll(i) + 0.5 * (s[i] + 1.0) * (s_macro_ur(i) - s_macro_ll(i));
267 }
268 Macro_elem_pt->macro_map(s_macro, x);
269 }
270
271 /// Global coordinates as function of local coordinates
272 /// at previous time "level" t (t=0: present; t>0: previous)
273 /// using the macro element representation
274 void get_x_from_macro_element(const unsigned& t,
275 const Vector<double>& s,
277 {
278 // Check that there is a macro element
279 if (Macro_elem_pt == 0)
280 {
281 throw OomphLibError("Macro Element pointer not set in this element\n",
282 OOMPH_CURRENT_FUNCTION,
283 OOMPH_EXCEPTION_LOCATION);
284 }
285
286 // Use the macro element representation
287 unsigned el_dim = dim();
288 Vector<double> s_macro(el_dim);
289 for (unsigned i = 0; i < el_dim; i++)
290 {
291 s_macro[i] =
292 s_macro_ll(i) + 0.5 * (s[i] + 1.0) * (s_macro_ur(i) - s_macro_ll(i));
293 }
294 Macro_elem_pt->macro_map(t, s_macro, x);
295 }
296
297 /// Return number of nodes on one face of the element. Always
298 /// nnode_1d^(el_dim - 1).
299 unsigned nnode_on_face() const
300 {
301 // c++ doesn't have pow(int, int) so we have to use all these casts...
302 return static_cast<unsigned>(std::pow(static_cast<double>(nnode_1d()),
303 static_cast<double>(dim() - 1)));
304 }
305
306 /// It's a Q element!
308 {
309 return ElementGeometry::Q;
310 }
311
312 private:
313 /// Pointer to vector of lower left vertex coords. in macro element
315
316 /// Pointer to vector of upper right vertex coords. in macro element
318 };
319
320
321 /// ///////////////////////////////////////////////////////////////////////
322 /// ///////////////////////////////////////////////////////////////////////
323 /// ///////////////////////////////////////////////////////////////////////
324
325
326 //========================================================================
327 /// Base class for Solid Qelements
328 //========================================================================
329 class QSolidElementBase : public virtual QElementBase,
330 public virtual SolidFiniteElement
331 {
332 public:
333 /// Constructor: Empty
335
336 /// Broken copy constructor
338
339 /// Broken assignment operator
340 /*void operator=(const QSolidElementBase&) = delete;*/
341
342 /// Set pointer to MacroElement -- overloads generic version
343 /// in RefineableQElement<2> and uses the MacroElement
344 /// also as the default for the "undeformed" configuration.
345 /// This assignment can/must be overwritten with
346 /// set_undeformed_macro_elem_pt(...) if the deformation of
347 /// the solid body is driven by a deformation of the
348 /// "current" Domain/MacroElement representation of it's boundary.
350 {
351 // Call the general Q version which sets up the storage
352 // for the reference coordinates
354 // Store pointer to macro element that represents the exact
355 // undeformed geomtry
357 }
358
359 /// Set pointers to "current" and "undeformed" MacroElements.
362 {
363 // Call the general Q version which sets up the storage
364 // for the reference coordinates
366 // Store pointer to macro element that represents the exact
367 // undeformed geomtry
369 }
370
371 /// Eulerian and Lagrangian coordinates as function of the
372 /// local coordinates: The Eulerian position is returned in
373 /// FE-interpolated form (\c x_fe) and then in the form obtained
374 /// from the "current" MacroElement representation (if it exists -- if not,
375 /// \c x is the same as \c x_fe). This allows the Domain/MacroElement-
376 /// based representation to be used to apply displacement boundary
377 /// conditions exactly. Ditto for the Lagrangian coordinates returned
378 /// in xi_fe and xi.
380 Vector<double>& x_fe,
382 Vector<double>& xi_fe,
383 Vector<double>& xi) const
384 {
385 // Lagrangian coordinate: Directly from
386 // underlying FE representation
387 unsigned n_xi = xi_fe.size();
388 for (unsigned i = 0; i < n_xi; i++)
389 {
390 xi_fe[i] = interpolated_xi(s, i);
391 }
392
393 // Lagrangian coordinate from FE representation again
395 {
396 unsigned n_xi = xi.size();
397 for (unsigned i = 0; i < n_xi; i++)
398 {
399 xi[i] = xi_fe[i];
400 }
401 }
402 // ...or refer to the "undeformed" MacroElement if it exists.
403 else
404 {
405 unsigned el_dim = dim();
406 Vector<double> s_macro(el_dim);
407 for (unsigned i = 0; i < el_dim; i++)
408 {
409 s_macro[i] = s_macro_ll(i) +
410 0.5 * (s[i] + 1.0) * (s_macro_ur(i) - s_macro_ll(i));
411 }
413 }
414
415
416 // Eulerian coordinate directly from underlying FE representation
417 unsigned n_x = x_fe.size();
418 for (unsigned i = 0; i < n_x; i++)
419 {
420 x_fe[i] = interpolated_x(s, i);
421 }
422
423 // Eulerian coordinate from FE representation again:
424 if (Macro_elem_pt == 0)
425 {
426 for (unsigned i = 0; i < n_x; i++)
427 {
428 x[i] = x_fe[i];
429 }
430 }
431 // or refer to the "current" MacroElement if it exists.
432 else
433 {
434 unsigned el_dim = dim();
435 Vector<double> s_macro(el_dim);
436 for (unsigned i = 0; i < el_dim; i++)
437 {
438 s_macro[i] = s_macro_ll(i) +
439 0.5 * (s[i] + 1.0) * (s_macro_ur(i) - s_macro_ll(i));
440 }
441 Macro_elem_pt->macro_map(s_macro, x);
442 }
443 }
444 };
445
446
447 /// ///////////////////////////////////////////////////////////////////////
448 /// ///////////////////////////////////////////////////////////////////////
449 /// ///////////////////////////////////////////////////////////////////////
450
451
452 //=======================================================================
453 /// General QElement class
454 ///
455 /// Empty, just establishes the template parameters
456 //=======================================================================
457 template<unsigned DIM, unsigned NNODE_1D>
459 {
460 };
461
462 //=======================================================================
463 /// Base class for all line elements
464 //=======================================================================
465 class LineElementBase : public virtual QElementBase
466 {
467 public:
468 /// Constructor. Empty
470
471 /// Number of vertex nodes in the element
472 virtual unsigned nvertex_node() const = 0;
473
474 /// Pointer to the j-th vertex node in the element
475 virtual Node* vertex_node_pt(const unsigned& j) const = 0;
476 };
477
478 //=======================================================================
479 /// General QElement class specialised to one spatial dimension
480 //=======================================================================
481 template<unsigned NNODE_1D>
482 class QElement<1, NNODE_1D> : public virtual LineElementBase
483 {
484 private:
485 /// Default integration rule: Gaussian integration of same 'order'
486 /// as the element
487 // This is sort of optimal, because it means that the integration is exact
488 // for the shape functions. Can overwrite this in specific element
489 // defintion.
491
492 public:
493 /// Constructor
495 {
496 // There are NNODE_1D nodes in this element
497 this->set_n_node(NNODE_1D);
498 // Set the dimensions of the element and the nodes, by default, both 1D
499 this->set_dimension(1);
500 // Assign pointer to default (full) integration_scheme
501 this->set_integration_scheme(&Default_integration_scheme);
502 }
503
504 /// Broken copy constructor
505 QElement(const QElement&) = delete;
506
507 /// Broken assignment operator
508 /*void operator=(const QElement&) = delete;*/
509
510 /// Calculate the geometric shape functions at local coordinate s
511 void shape(const Vector<double>& s, Shape& psi) const;
512
513 /// Compute the geometric shape functions and
514 /// derivatives w.r.t. local coordinates at local coordinate s
515 void dshape_local(const Vector<double>& s,
516 Shape& psi,
517 DShape& dpsids) const;
518
519 /// Compute the geometric shape functions, derivatives and
520 /// second derivatives w.r.t. local coordinates at local coordinate s
521 /// d2psids(i,0) = \f$ d^2 \psi_j / d s^2 \f$
522 void d2shape_local(const Vector<double>& s,
523 Shape& psi,
524 DShape& dpsids,
525 DShape& d2psids) const;
526
527 /// Overload the template-free interface for the calculation of
528 /// inverse jacobian matrix. This is a one-dimensional element, so
529 /// use the 1D version.
531 DenseMatrix<double>& inverse_jacobian) const
532 {
533 return FiniteElement::invert_jacobian<1>(jacobian, inverse_jacobian);
534 }
535
536 /// Min. value of local coordinate
537 double s_min() const
538 {
539 return -1.0;
540 }
541
542 /// Max. value of local coordinate
543 double s_max() const
544 {
545 return 1.0;
546 }
547
548 /// Number of vertex nodes in the element
549 unsigned nvertex_node() const
550 {
551 return 2;
552 }
553
554 /// Pointer to the j-th vertex node in the element
555 Node* vertex_node_pt(const unsigned& j) const
556 {
557 unsigned n_node_1d = nnode_1d();
558 Node* nod_pt;
559 switch (j)
560 {
561 case 0:
562 nod_pt = node_pt(0);
563 break;
564 case 1:
565 nod_pt = node_pt(n_node_1d - 1);
566 break;
567 default:
568 std::ostringstream error_message;
569 error_message << "Vertex node number is " << j
570 << " but must be from 0 to 1\n";
571
572 throw OomphLibError(error_message.str(),
573 OOMPH_CURRENT_FUNCTION,
574 OOMPH_EXCEPTION_LOCATION);
575 }
576 return nod_pt;
577 }
578
579 /// Get local coordinates of node j in the element; vector sets its own size
580 void local_coordinate_of_node(const unsigned& j, Vector<double>& s) const
581 {
582 s.resize(1);
583 s[0] = this->s_min() +
584 double(j) / double(NNODE_1D - 1) * (this->s_max() - this->s_min());
585 }
586
587
588 /// Get the local fraction of node j in the element
589 void local_fraction_of_node(const unsigned& j, Vector<double>& s_fraction)
590 {
591 s_fraction.resize(1);
592 s_fraction[0] = double(j) / double(NNODE_1D - 1);
593 }
594
595
596 /// This function returns the local fraction of all nodes at the n-th
597 /// position in a one dimensional expansion along the i-th local coordinate
598 inline double local_one_d_fraction_of_node(const unsigned& n1d,
599 const unsigned& i)
600 {
601 // It's just the value of the node divided by the number of 1-D nodes
602 return double(n1d) / double(NNODE_1D - 1);
603 }
604
605 /// Get the node at the specified local coordinate
606 Node* get_node_at_local_coordinate(const Vector<double>& s) const;
607
608 /// Number of nodes along each element edge
609 unsigned nnode_1d() const
610 {
611 return NNODE_1D;
612 }
613
614 /// Return the number of actual plot points for paraview
615 /// plot with parameter nplot.
616 unsigned nplot_points_paraview(const unsigned& nplot) const
617 {
618 return nplot;
619 }
620
621 /// Return the number of local sub-elements for paraview plot with
622 /// parameter nplot.
623 unsigned nsub_elements_paraview(const unsigned& nplot) const
624 {
625 return (nplot - 1);
626 }
627
628 /// Fill in the offset information for paraview plot.
629 /// Needs to be implemented for each new geometric element type; see
630 /// http://www.vtk.org/VTK/img/file-formats.pdf
631 void write_paraview_output_offset_information(std::ofstream& file_out,
632 const unsigned& nplot,
633 unsigned& counter) const
634 {
635 // Number of local elements we want to plot over
636 unsigned plot = nsub_elements_paraview(nplot);
637
638 // loops over the i-th local element in parent element
639 for (unsigned i = 0; i < plot; i++)
640 {
641 file_out << i + counter << " " << i + 1 + counter << std::endl;
642 }
643 counter += nplot_points_paraview(nplot);
644 }
645
646 /// Return the paraview element type.
647 /// Needs to be implemented for each new geometric element type; see
648 /// http://www.vtk.org/VTK/img/file-formats.pdf
649 /// Use type "VTK_LINE" (== 3) for 2D quad elements
650 void write_paraview_type(std::ofstream& file_out,
651 const unsigned& nplot) const
652 {
653 unsigned local_loop = nsub_elements_paraview(nplot);
654 for (unsigned i = 0; i < local_loop; i++)
655 {
656 file_out << "3" << std::endl;
657 }
658 }
659
660 /// Return the offsets for the paraview sub-elements. Needs
661 /// to be implemented for each new geometric element type; see
662 /// http://www.vtk.org/VTK/img/file-formats.pdf
663 void write_paraview_offsets(std::ofstream& file_out,
664 const unsigned& nplot,
665 unsigned& offset_sum) const
666 {
667 // Loop over all local elements and add its offset to the overall
668 // offset_sum
669 unsigned local_loop = nsub_elements_paraview(nplot);
670 for (unsigned i = 0; i < local_loop; i++)
671 {
672 offset_sum += 2;
673 file_out << offset_sum << std::endl;
674 }
675 }
676
677 /// Output
678 void output(std::ostream& outfile);
679
680 /// Output at n_plot points
681 void output(std::ostream& outfile, const unsigned& n_plot);
682
683 /// C-style output
684 void output(FILE* file_pt);
685
686 /// C_style output at n_plot points
687 void output(FILE* file_pt, const unsigned& n_plot);
688
689
690 /// Get cector of local coordinates of plot point i (when plotting
691 /// nplot points in each "coordinate direction).
693 const unsigned& i,
694 const unsigned& nplot,
696 const bool& use_equally_spaced_interior_sample_points = false) const
697 {
698 if (nplot > 1)
699 {
700 s[0] = -1.0 + 2.0 * double(i) / double(nplot - 1);
701 if (use_equally_spaced_interior_sample_points)
702 {
703 double range = 2.0;
704 double dx_new = range / double(nplot);
705 double range_new = double(nplot - 1) * dx_new;
706 s[0] = -1.0 + 0.5 * dx_new + range_new * (1.0 + s[0]) / range;
707 }
708 }
709 else
710 {
711 s[0] = 0.0;
712 }
713 }
714
715 /// Return string for tecplot zone header (when plotting
716 /// nplot points in each "coordinate direction)
717 std::string tecplot_zone_string(const unsigned& nplot) const
718 {
719 std::ostringstream header;
720 header << "ZONE I=" << nplot << "\n";
721 return header.str();
722 }
723
724 /// Return total number of plot points (when plotting
725 /// nplot points in each "coordinate direction)
726 unsigned nplot_points(const unsigned& nplot) const
727 {
728 unsigned DIM = 1;
729 unsigned np = 1;
730 for (unsigned i = 0; i < DIM; i++)
731 {
732 np *= nplot;
733 }
734 return np;
735 }
736
737 /// Get the number of the ith node on face face_index in the bulk node
738 /// vector.
739 unsigned get_bulk_node_number(const int& face_index,
740 const unsigned& i) const
741 {
742 face_node_number_error_check(i);
743
744 if (face_index == -1)
745 {
746 return 0;
747 }
748 else if (face_index == +1)
749 {
750 return nnode_1d() - 1;
751 }
752 else
753 {
754 std::string err = "Face index should be in {-1, +1}.";
755 throw OomphLibError(
756 err, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
757 }
758 }
759
760 /// Get the sign of the outer unit normal on the face given by face_index.
761 int face_outer_unit_normal_sign(const int& face_index) const
762 {
763#ifdef PARANOID
764 if (std::abs(face_index) != 1)
765 {
766 std::string err = "Face index should be in {-1, +1}.";
767 throw OomphLibError(
768 err, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
769 }
770#endif
771 return face_index;
772 }
773
774 /// Get a pointer to the function mapping face coordinates to bulk
775 /// coordinates
777 const int& face_index) const
778 {
779 if (face_index == 1)
780 {
782 }
783 else if (face_index == -1)
784 {
786 }
787 else
788 {
789 std::string err = "Face index should be in {-1, +1}.";
790 throw OomphLibError(
791 err, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
792 }
793 }
794
795 /// Get a pointer to the derivative of the mapping from face to bulk
796 /// coordinates.
798 const int& face_index) const
799 {
800 if (face_index == 1)
801 {
803 }
804 else if (face_index == -1)
805 {
807 }
808 else
809 {
810 std::string err = "Face index should be in {-1, +1}.";
811 throw OomphLibError(
812 err, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
813 }
814 }
815 };
816
817 //=======================================================================
818 /// Base class for all quad elements
819 //=======================================================================
820 class QuadElementBase : public virtual QElementBase
821 {
822 public:
823 /// Constructor. Empty
825
826 /// Number of vertex nodes in the element
827 virtual unsigned nvertex_node() const = 0;
828
829 /// Pointer to the j-th vertex node in the element
830 virtual Node* vertex_node_pt(const unsigned& j) const = 0;
831 };
832
833 //=======================================================================
834 /// General QElement class specialised to two spatial dimensions
835 //=======================================================================
836 template<unsigned NNODE_1D>
837 class QElement<2, NNODE_1D> : public virtual QuadElementBase
838 {
839 private:
840 /// Default integration rule: Gaussian integration of same 'order'
841 /// as the element
842 // N.B. This is sort of optimal, because it means that the integration is
843 // exact for the shape functions. Can overwrite this in specific element
844 // defintion
846
847 public:
848 /// Constructor
850 {
851 // There are NNODE_1D*NNODE_1D nodes in this element
852 this->set_n_node(NNODE_1D * NNODE_1D);
853 // Set the dimensions of the element and the nodes, by default, both 2D
854 set_dimension(2);
855 // Assign default (full) spatial integration scheme
856 set_integration_scheme(&Default_integration_scheme);
857 }
858
859 /// Broken copy constructor
860 QElement(const QElement&) = delete;
861
862 /// Broken assignment operator
863 /*void operator=(const QElement&) = delete;*/
864
865 /// Calculate the geometric shape functions at local coordinate s
866 void shape(const Vector<double>& s, Shape& psi) const;
867
868 /// Compute the geometric shape functions and
869 /// derivatives w.r.t. local coordinates at local coordinate s
870 void dshape_local(const Vector<double>& s,
871 Shape& psi,
872 DShape& dpsids) const;
873
874 /// Compute the geometric shape functions, derivatives and
875 /// second derivatives w.r.t. local coordinates at local coordinate s
876 /// d2psids(i,0) = \f$ \partial^2 \psi_j / \partial s_0^2 \f$
877 /// d2psids(i,1) = \f$ \partial^2 \psi_j / \partial s_1^2 \f$
878 /// d2psids(i,2) = \f$ \partial^2 \psi_j / \partial s_0 \partial s_1 \f$
879 void d2shape_local(const Vector<double>& s,
880 Shape& psi,
881 DShape& dpsids,
882 DShape& d2psids) const;
883
884 /// Overload the template-free interface for the calculation of
885 /// inverse jacobian matrix. This is a two-dimensional element, so use
886 /// the two-d version.
888 DenseMatrix<double>& inverse_jacobian) const
889 {
890 return FiniteElement::invert_jacobian<2>(jacobian, inverse_jacobian);
891 }
892
893 /// Min. value of local coordinate
894 double s_min() const
895 {
896 return -1.0;
897 }
898
899 /// Max. value of local coordinate
900 double s_max() const
901 {
902 return 1.0;
903 }
904
905
906 /// Number of vertex nodes in the element
907 unsigned nvertex_node() const
908 {
909 return 4;
910 }
911
912 /// Pointer to the j-th vertex node in the element
913 Node* vertex_node_pt(const unsigned& j) const
914 {
915 unsigned n_node_1d = nnode_1d();
916 Node* nod_pt;
917 switch (j)
918 {
919 case 0:
920 nod_pt = node_pt(0);
921 break;
922 case 1:
923 nod_pt = node_pt(n_node_1d - 1);
924 break;
925 case 2:
926 nod_pt = node_pt(n_node_1d * (n_node_1d - 1));
927 break;
928 case 3:
929 nod_pt = node_pt(n_node_1d * n_node_1d - 1);
930 break;
931 default:
932 std::ostringstream error_message;
933 error_message << "Vertex node number is " << j
934 << " but must be from 0 to 3\n";
935
936 throw OomphLibError(error_message.str(),
937 OOMPH_CURRENT_FUNCTION,
938 OOMPH_EXCEPTION_LOCATION);
939 }
940 return nod_pt;
941 }
942
943
944 /// Get local coordinates of node j in the element; vector sets its own size
945 void local_coordinate_of_node(const unsigned& j, Vector<double>& s) const
946 {
947 s.resize(2);
948 unsigned j0 = j % NNODE_1D;
949 unsigned j1 = j / NNODE_1D;
950 const double S_min = this->s_min();
951 const double S_range = this->s_max() - S_min;
952 s[0] = S_min + double(j0) / double(NNODE_1D - 1) * S_range;
953 s[1] = S_min + double(j1) / double(NNODE_1D - 1) * S_range;
954 }
955
956 /// Get the local fraction of node j in the element
957 void local_fraction_of_node(const unsigned& j, Vector<double>& s_fraction)
958 {
959 s_fraction.resize(2);
960 unsigned j0 = j % NNODE_1D;
961 unsigned j1 = j / NNODE_1D;
962 s_fraction[0] = double(j0) / double(NNODE_1D - 1);
963 s_fraction[1] = double(j1) / double(NNODE_1D - 1);
964 }
965
966 /// This function returns the local fraction of ant nodes in the n-th
967 /// positoin in a one dimensional expansion along the i-th local coordinate
968 inline double local_one_d_fraction_of_node(const unsigned& n1d,
969 const unsigned& i)
970 {
971 // It's just the value of the node divided by the number of 1-D nodes
972 return double(n1d) / double(NNODE_1D - 1);
973 }
974
975 /// Get the node at the specified local coordinate
976 Node* get_node_at_local_coordinate(const Vector<double>& s) const;
977
978 /// Number of nodes along each element edge
979 unsigned nnode_1d() const
980 {
981 return NNODE_1D;
982 }
983
984 /// Return the number of actual plot points for paraview
985 /// plot with parameter nplot.
986 unsigned nplot_points_paraview(const unsigned& nplot) const
987 {
988 return nplot * nplot;
989 }
990
991 /// Return the number of local sub-elements for paraview plot with
992 /// parameter nplot.
993 unsigned nsub_elements_paraview(const unsigned& nplot) const
994 {
995 return (nplot - 1) * (nplot - 1);
996 }
997
998 /// Fill in the offset information for paraview plot.
999 /// Needs to be implemented for each new geometric element type; see
1000 /// http://www.vtk.org/VTK/img/file-formats.pdf
1001 void write_paraview_output_offset_information(std::ofstream& file_out,
1002 const unsigned& nplot,
1003 unsigned& counter) const
1004 {
1005 // Number of local elements we want to plot over
1006 unsigned plot = nsub_elements_paraview(nplot);
1007
1008 // loops over the i-th local element in parent element
1009 for (unsigned i = 0; i < plot; i++)
1010 {
1011 unsigned d = (i - (i % (nplot - 1))) / (nplot - 1);
1012
1013 file_out << i % (nplot - 1) + d * nplot + counter << " "
1014 << i % (nplot - 1) + 1 + d * nplot + counter << " "
1015 << i % (nplot - 1) + 1 + (d + 1) * nplot + counter << " "
1016 << i % (nplot - 1) + (d + 1) * nplot + counter << std::endl;
1017 }
1018 counter += nplot_points_paraview(nplot);
1019 }
1020
1021 /// Return the paraview element type.
1022 /// Needs to be implemented for each new geometric element type; see
1023 /// http://www.vtk.org/VTK/img/file-formats.pdf
1024 /// Use type "VTK_QUAD" (== 9) for 2D quad elements
1025 void write_paraview_type(std::ofstream& file_out,
1026 const unsigned& nplot) const
1027 {
1028 unsigned local_loop = nsub_elements_paraview(nplot);
1029 for (unsigned i = 0; i < local_loop; i++)
1030 {
1031 file_out << "9" << std::endl;
1032 }
1033 }
1034
1035 /// Return the offsets for the paraview sub-elements. Needs
1036 /// to be implemented for each new geometric element type; see
1037 /// http://www.vtk.org/VTK/img/file-formats.pdf
1038 void write_paraview_offsets(std::ofstream& file_out,
1039 const unsigned& nplot,
1040 unsigned& offset_sum) const
1041 {
1042 // Loop over all local elements and add its offset to the overall
1043 // offset_sum
1044 unsigned local_loop = nsub_elements_paraview(nplot);
1045 for (unsigned i = 0; i < local_loop; i++)
1046 {
1047 offset_sum += 4;
1048 file_out << offset_sum << std::endl;
1049 }
1050 }
1051
1052 /// Output
1053 void output(std::ostream& outfile);
1054
1055 /// Output at n_plot points
1056 void output(std::ostream& outfile, const unsigned& n_plot);
1057
1058 /// C-style output
1059 void output(FILE* file_pt);
1060
1061 /// C_style output at n_plot points
1062 void output(FILE* file_pt, const unsigned& n_plot);
1063
1064
1065 /// Get cector of local coordinates of plot point i (when plotting
1066 /// nplot points in each "coordinate direction).
1068 const unsigned& i,
1069 const unsigned& nplot,
1071 const bool& use_equally_spaced_interior_sample_points = false) const
1072 {
1073 if (nplot > 1)
1074 {
1075 unsigned i0 = i % nplot;
1076 unsigned i1 = (i - i0) / nplot;
1077
1078 s[0] = -1.0 + 2.0 * double(i0) / double(nplot - 1);
1079 s[1] = -1.0 + 2.0 * double(i1) / double(nplot - 1);
1080 if (use_equally_spaced_interior_sample_points)
1081 {
1082 double range = 2.0;
1083 double dx_new = range / double(nplot);
1084 double range_new = double(nplot - 1) * dx_new;
1085 s[0] = -1.0 + 0.5 * dx_new + range_new * (1.0 + s[0]) / range;
1086 s[1] = -1.0 + 0.5 * dx_new + range_new * (1.0 + s[1]) / range;
1087 }
1088 }
1089 else
1090 {
1091 s[0] = 0.0;
1092 s[1] = 0.0;
1093 }
1094 }
1095
1096 /// Return string for tecplot zone header (when plotting
1097 /// nplot points in each "coordinate direction)
1098 std::string tecplot_zone_string(const unsigned& nplot) const
1099 {
1100 std::ostringstream header;
1101 header << "ZONE I=" << nplot << ", J=" << nplot << "\n";
1102 return header.str();
1103 }
1104
1105 /// Return total number of plot points (when plotting
1106 /// nplot points in each "coordinate direction)
1107 unsigned nplot_points(const unsigned& nplot) const
1108 {
1109 unsigned DIM = 2;
1110 unsigned np = 1;
1111 for (unsigned i = 0; i < DIM; i++)
1112 {
1113 np *= nplot;
1114 }
1115 return np;
1116 };
1117
1118 /// Get the number of the ith node on face face_index (in the bulk node
1119 /// vector).
1120 unsigned get_bulk_node_number(const int& face_index,
1121 const unsigned& i) const
1122 {
1123 face_node_number_error_check(i);
1124
1125 const unsigned nn1d = nnode_1d();
1126
1127 if (face_index == -1)
1128 {
1129 return i * nn1d;
1130 }
1131 else if (face_index == +1)
1132 {
1133 return nn1d * i + nn1d - 1;
1134 }
1135 else if (face_index == -2)
1136 {
1137 return i;
1138 }
1139 else if (face_index == +2)
1140 {
1141 return nn1d * (nn1d - 1) + i;
1142 }
1143 else
1144 {
1145 std::string err = "Face index should be in {-1, +1, -2, +2}.";
1146 throw OomphLibError(
1147 err, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
1148 }
1149 }
1150
1151 /// Get the sign of the outer unit normal on the face given by face_index.
1152 int face_outer_unit_normal_sign(const int& face_index) const
1153 {
1154 if (face_index < 0)
1155 {
1156 return 1;
1157 }
1158 else if (face_index > 0)
1159 {
1160 return -1;
1161 }
1162 else
1163 {
1164 std::string err = "Face index should be one of {-1, +1, -2, +2}.";
1165 throw OomphLibError(
1166 err, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
1167 }
1168 }
1169
1170 /// Get a pointer to the function mapping face coordinates to bulk
1171 /// coordinates
1173 const int& face_index) const
1174 {
1175 if (face_index == 1)
1176 {
1178 }
1179 else if (face_index == -1)
1180 {
1182 }
1183 else if (face_index == -2)
1184 {
1186 }
1187 else if (face_index == 2)
1188 {
1190 }
1191 else
1192 {
1193 std::string err = "Face index should be in {-1, +1}.";
1194 throw OomphLibError(
1195 err, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
1196 }
1197 }
1198
1199 /// Get a pointer to the derivative of the mapping from face to bulk
1200 /// coordinates.
1202 const int& face_index) const
1203 {
1204 if (face_index == 1)
1205 {
1207 }
1208 else if (face_index == -1)
1209 {
1211 }
1212 else if (face_index == -2)
1213 {
1215 }
1216 else if (face_index == 2)
1217 {
1219 }
1220 else
1221 {
1222 std::string err = "Face index should be in {-1, +1}.";
1223 throw OomphLibError(
1224 err, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
1225 }
1226 }
1227 };
1228
1229 //=======================================================================
1230 /// Base class for all brick elements
1231 //=======================================================================
1232 class BrickElementBase : public virtual QElementBase
1233 {
1234 public:
1235 /// Constructor. Empty
1237
1238 /// Number of vertex nodes in the element
1239 virtual unsigned nvertex_node() const = 0;
1240
1241 /// Pointer to the j-th vertex node in the element
1242 virtual Node* vertex_node_pt(const unsigned& j) const = 0;
1243 };
1244
1245 //=======================================================================
1246 /// General QElement class specialised to three spatial dimensions
1247 //=======================================================================
1248 template<unsigned NNODE_1D>
1249 class QElement<3, NNODE_1D> : public virtual BrickElementBase
1250 {
1251 private:
1252 /// Default integration rule: Gaussian integration of same 'order'
1253 /// as the element
1254 // N.B. This is sort of optimal, because it means that the integration is
1255 // exact for the shape functions. Can overwrite this in specific element
1256 // defintion
1258
1259 public:
1260 /// Constructor
1262 {
1263 // There are NNODE_1D^3 nodes in this element
1264 this->set_n_node(NNODE_1D * NNODE_1D * NNODE_1D);
1265 // Set the dimensions of the element and the nodes, by default, both 3D
1266 set_dimension(3);
1267 // Assign default (full_ spatial integration_scheme
1268 set_integration_scheme(&Default_integration_scheme);
1269 }
1270
1271
1272 /// Broken copy constructor
1273 QElement(const QElement&) = delete;
1274
1275 /// Broken assignment operator
1276 /*void operator=(const QElement&) = delete;*/
1277
1278 /// Calculate the geometric shape functions at local coordinate s
1279 void shape(const Vector<double>& s, Shape& psi) const;
1280
1281 /// Compute the geometric shape functions and
1282 /// derivatives w.r.t. local coordinates at local coordinate s
1283 void dshape_local(const Vector<double>& s,
1284 Shape& psi,
1285 DShape& dpsids) const;
1286
1287 /// Compute the geometric shape functions, derivatives and
1288 /// second derivatives w.r.t. local coordinates at local coordinate s.
1289 /// d2psids(i,0) = \f$ \partial^2 \psi_j / \partial s_0^2 \f$
1290 /// d2psids(i,1) = \f$ \partial^2 \psi_j / \partial s_1^2 \f$
1291 /// d2psids(i,2) = \f$ \partial^2 \psi_j / \partial s_2^2 \f$
1292 /// d2psids(i,3) = \f$ \partial^2 \psi_j / \partial s_0 \partial s_1 \f$
1293 /// d2psids(i,4) = \f$ \partial^2 \psi_j / \partial s_0 \partial s_2 \f$
1294 /// d2psids(i,5) = \f$ \partial^2 \psi_j / \partial s_1 \partial s_2 \f$
1295 void d2shape_local(const Vector<double>& s,
1296 Shape& psi,
1297 DShape& dpsids,
1298 DShape& d2psids) const;
1299
1300
1301 /// Overload the template-free interface for the calculation of
1302 /// the inverse jacobian mapping. This is a three-dimensional element,
1303 /// so use the 3d version
1305 DenseMatrix<double>& inverse_jacobian) const
1306 {
1307 return FiniteElement::invert_jacobian<3>(jacobian, inverse_jacobian);
1308 }
1309
1310 /// Min. value of local coordinate
1311 double s_min() const
1312 {
1313 return -1.0;
1314 }
1315
1316 /// Max. value of local coordinate
1317 double s_max() const
1318 {
1319 return 1.0;
1320 }
1321
1322 /// Number of vertex nodes in the element
1323 unsigned nvertex_node() const
1324 {
1325 return 8;
1326 }
1327
1328 /// Pointer to the j-th vertex node in the element
1329 Node* vertex_node_pt(const unsigned& j) const
1330 {
1331 unsigned N = nnode_1d();
1332 Node* nod_pt;
1333 switch (j)
1334 {
1335 case 0:
1336 nod_pt = node_pt(0);
1337 break;
1338 case 1:
1339 nod_pt = node_pt(N - 1);
1340 break;
1341 case 2:
1342 nod_pt = node_pt(N * (N - 1));
1343 break;
1344 case 3:
1345 nod_pt = node_pt(N * N - 1);
1346 break;
1347 case 4:
1348 nod_pt = node_pt(N * N * (N - 1));
1349 break;
1350 case 5:
1351 nod_pt = node_pt(N * N * (N - 1) + (N - 1));
1352 break;
1353 case 6:
1354 nod_pt = node_pt(N * N * N - N);
1355 break;
1356 case 7:
1357 nod_pt = node_pt(N * N * N - 1);
1358 break;
1359 default:
1360 std::ostringstream error_message;
1361 error_message << "Vertex node number is " << j
1362 << " but must be from 0 to 7\n";
1363
1364 throw OomphLibError(error_message.str(),
1365 OOMPH_CURRENT_FUNCTION,
1366 OOMPH_EXCEPTION_LOCATION);
1367 }
1368 return nod_pt;
1369 }
1370
1371 /// Get local coordinates of node j in the element; vector sets its own size
1372 void local_coordinate_of_node(const unsigned& j, Vector<double>& s) const
1373 {
1374 s.resize(3);
1375 unsigned j0 = j % NNODE_1D;
1376 unsigned j1 = (j / NNODE_1D) % NNODE_1D;
1377 unsigned j2 = j / (NNODE_1D * NNODE_1D);
1378 const double S_min = this->s_min();
1379 const double S_range = this->s_max() - S_min;
1380
1381 s[0] = S_min + double(j0) / double(NNODE_1D - 1) * S_range;
1382 s[1] = S_min + double(j1) / double(NNODE_1D - 1) * S_range;
1383 s[2] = S_min + double(j2) / double(NNODE_1D - 1) * S_range;
1384 }
1385
1386 /// Get the local fraction of node j in the element
1387 void local_fraction_of_node(const unsigned& j, Vector<double>& s_fraction)
1388 {
1389 s_fraction.resize(3);
1390 unsigned j0 = j % NNODE_1D;
1391 unsigned j1 = (j / NNODE_1D) % NNODE_1D;
1392 unsigned j2 = j / (NNODE_1D * NNODE_1D);
1393 s_fraction[0] = double(j0) / double(NNODE_1D - 1);
1394 s_fraction[1] = double(j1) / double(NNODE_1D - 1);
1395 s_fraction[2] = double(j2) / double(NNODE_1D - 1);
1396 }
1397
1398 /// This function returns the local fraction of any nodes in the n-th
1399 /// positoin in a one dimensional expansion along the i-th local coordinate
1400 inline double local_one_d_fraction_of_node(const unsigned& n1d,
1401 const unsigned& i)
1402 {
1403 // It's just the value of the node divided by the number of 1-D nodes
1404 return double(n1d) / double(NNODE_1D - 1);
1405 }
1406
1407 /// Get the node at the specified local coordinate
1408 Node* get_node_at_local_coordinate(const Vector<double>& s) const;
1409
1410 /// Number of nodes along each element edge
1411 unsigned nnode_1d() const
1412 {
1413 return NNODE_1D;
1414 }
1415
1416 /// Return the number of actual plot points for paraview
1417 /// plot with parameter nplot.
1418 unsigned nplot_points_paraview(const unsigned& nplot) const
1419 {
1420 return nplot * nplot * nplot;
1421 }
1422
1423 /// Return the number of local sub-elements for paraview plot with
1424 /// parameter nplot.
1425 unsigned nsub_elements_paraview(const unsigned& nplot) const
1426 {
1427 return (nplot - 1) * (nplot - 1) * (nplot - 1);
1428 }
1429
1430 /// Fill in the offset information for paraview plot.
1431 /// Needs to be implemented for each new geometric element type; see
1432 /// http://www.vtk.org/VTK/img/file-formats.pdf
1433 void write_paraview_output_offset_information(std::ofstream& file_out,
1434 const unsigned& nplot,
1435 unsigned& counter) const
1436 {
1437 // Number of local elements we want to plot over
1438 unsigned plot = nsub_elements_paraview(nplot);
1439
1440 for (unsigned j = 0; j < plot; j += (nplot - 1) * (nplot - 1) + 1)
1441 {
1442 // To keep track of how many cross-sections we've looped over
1443 unsigned r = ((j - (j % ((nplot - 1) * (nplot - 1)))) /
1444 ((nplot - 1) * (nplot - 1)));
1445
1446 // loop over all the elemnets in this sublevel
1447 unsigned sub_plot = (nplot - 1) * (nplot - 1);
1448
1449 // loops over the i-th local element in parent element
1450 for (unsigned i = 0; i < sub_plot; i++)
1451 {
1452 unsigned d = ((i - (i % (nplot - 1))) / (nplot - 1));
1453
1454
1455 // Lower level of rectangle
1456 file_out
1457 << i % (nplot - 1) + d * nplot + r * nplot * nplot + counter << " "
1458 << i % (nplot - 1) + 1 + d * nplot + r * nplot * nplot + counter
1459 << " "
1460 << i % (nplot - 1) + 1 + (d + 1) * nplot + r * nplot * nplot +
1461 counter
1462 << " "
1463 << i % (nplot - 1) + (d + 1) * nplot + r * nplot * nplot + counter
1464 << " "
1465
1466 // Upper level of rectangle
1467 << i % (nplot - 1) + d * nplot + (r + 1) * nplot * nplot + counter
1468 << " "
1469 << i % (nplot - 1) + 1 + d * nplot + (r + 1) * nplot * nplot +
1470 counter
1471 << " "
1472 << i % (nplot - 1) + 1 + (d + 1) * nplot + (r + 1) * nplot * nplot +
1473 counter
1474 << " "
1475 << i % (nplot - 1) + (d + 1) * nplot + (r + 1) * nplot * nplot +
1476 counter
1477 << std::endl;
1478 }
1479 }
1480 counter += nplot_points_paraview(nplot);
1481 }
1482
1483 /// Return the paraview element type.
1484 /// Needs to be implemented for each new geometric element type; see
1485 /// http://www.vtk.org/VTK/img/file-formats.pdf
1486 /// Use type "VTK_HEXAHEDRON" (== 12) for 2D quad elements
1487 void write_paraview_type(std::ofstream& file_out,
1488 const unsigned& nplot) const
1489 {
1490 unsigned local_loop = nsub_elements_paraview(nplot);
1491 for (unsigned i = 0; i < local_loop; i++)
1492 {
1493 file_out << "12" << std::endl;
1494 }
1495 }
1496
1497 /// Return the offsets for the paraview sub-elements. Needs
1498 /// to be implemented for each new geometric element type; see
1499 /// http://www.vtk.org/VTK/img/file-formats.pdf
1500 void write_paraview_offsets(std::ofstream& file_out,
1501 const unsigned& nplot,
1502 unsigned& offset_sum) const
1503 {
1504 unsigned local_loop = nsub_elements_paraview(nplot);
1505 for (unsigned i = 0; i < local_loop; i++)
1506 {
1507 offset_sum += 8;
1508 file_out << offset_sum << std::endl;
1509 }
1510 }
1511
1512 /// Output
1513 void output(std::ostream& outfile);
1514
1515 /// Output at n_plot points
1516 void output(std::ostream& outfile, const unsigned& n_plot);
1517
1518 /// C-style output
1519 void output(FILE* file_pt);
1520
1521 /// C_style output at n_plot points
1522 void output(FILE* file_pt, const unsigned& n_plot);
1523
1524 /// Get cector of local coordinates of plot point i (when plotting
1525 /// nplot points in each "coordinate direction).
1527 const unsigned& i,
1528 const unsigned& nplot,
1530 const bool& use_equally_spaced_interior_sample_points = false) const
1531 {
1532 if (nplot > 1)
1533 {
1534 unsigned i01 = i % (nplot * nplot);
1535 unsigned i0 = i01 % nplot;
1536 unsigned i1 = (i01 - i0) / nplot;
1537 unsigned i2 = (i - i01) / (nplot * nplot);
1538
1539 s[0] = -1.0 + 2.0 * double(i0) / double(nplot - 1);
1540 s[1] = -1.0 + 2.0 * double(i1) / double(nplot - 1);
1541 s[2] = -1.0 + 2.0 * double(i2) / double(nplot - 1);
1542 if (use_equally_spaced_interior_sample_points)
1543 {
1544 double range = 2.0;
1545 double dx_new = range / double(nplot);
1546 double range_new = double(nplot - 1) * dx_new;
1547 s[0] = -1.0 + 0.5 * dx_new + range_new * (1.0 + s[0]) / range;
1548 s[1] = -1.0 + 0.5 * dx_new + range_new * (1.0 + s[1]) / range;
1549 s[2] = -1.0 + 0.5 * dx_new + range_new * (1.0 + s[2]) / range;
1550 }
1551 }
1552 else
1553 {
1554 s[0] = 0.0;
1555 s[1] = 0.0;
1556 s[2] = 0.0;
1557 }
1558 }
1559
1560 /// Return string for tecplot zone header (when plotting
1561 /// nplot points in each "coordinate direction)
1562 std::string tecplot_zone_string(const unsigned& nplot) const
1563 {
1564 std::ostringstream header;
1565 header << "ZONE I=" << nplot << ", J=" << nplot << ", K=" << nplot
1566 << "\n";
1567 return header.str();
1568 }
1569
1570 /// Return total number of plot points (when plotting
1571 /// nplot points in each "coordinate direction)
1572 unsigned nplot_points(const unsigned& nplot) const
1573 {
1574 unsigned DIM = 3;
1575 unsigned np = 1;
1576 for (unsigned i = 0; i < DIM; i++)
1577 {
1578 np *= nplot;
1579 }
1580 return np;
1581 };
1582
1583 /// Get the number of the ith node on face face_index in the bulk node
1584 /// vector.
1585 unsigned get_bulk_node_number(const int& face_index,
1586 const unsigned& i) const
1587 {
1588 face_node_number_error_check(i);
1589
1590 const unsigned nn1d = nnode_1d();
1591
1592 if (face_index == -1)
1593 {
1594 return i * nn1d;
1595 }
1596 else if (face_index == +1)
1597 {
1598 return i * nn1d + (nn1d - 1);
1599 }
1600 else if (face_index == -2)
1601 {
1602 return i % nn1d + (i / nn1d) * nn1d * nn1d;
1603 }
1604 else if (face_index == +2)
1605 {
1606 return i % nn1d + (i / nn1d) * nn1d * nn1d + (nn1d * (nn1d - 1));
1607 }
1608 else if (face_index == -3)
1609 {
1610 return i;
1611 }
1612 else if (face_index == +3)
1613 {
1614 return i + (nn1d * nn1d) * (nn1d - 1);
1615 }
1616 else
1617 {
1618 std::string err = "Face index should be in {-1, +1, -2, +2, -3, +3}.";
1619 throw OomphLibError(
1620 err, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
1621 }
1622 }
1623
1624 /// Get the sign of the outer unit normal on the face given by face_index.
1625 int face_outer_unit_normal_sign(const int& face_index) const
1626 {
1627 if (face_index == -3)
1628 {
1629 return -1;
1630 }
1631 else if (face_index == +3)
1632 {
1633 return 1;
1634 }
1635 else if (face_index == -2)
1636 {
1637 return 1;
1638 }
1639 else if (face_index == 2)
1640 {
1641 return -1;
1642 }
1643 else if (face_index == -1)
1644 {
1645 return -1;
1646 }
1647 else if (face_index == 1)
1648 {
1649 return 1;
1650 }
1651 else
1652 {
1653 std::string err = "Face index should be in {-1, +1, -2, +2, -3, +3}.";
1654 throw OomphLibError(
1655 err, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
1656 }
1657 }
1658
1659 /// Get a pointer to the function mapping face coordinates to bulk
1660 /// coordinates
1662 const int& face_index) const
1663 {
1664 if (face_index == 1)
1665 {
1667 }
1668 else if (face_index == -1)
1669 {
1671 }
1672 else if (face_index == -2)
1673 {
1675 }
1676 else if (face_index == 2)
1677 {
1679 }
1680 else if (face_index == -3)
1681 {
1683 }
1684 else if (face_index == 3)
1685 {
1687 }
1688 else
1689 {
1690 std::string err = "Face index should be in {-1, +1}.";
1691 throw OomphLibError(
1692 err, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
1693 }
1694 }
1695
1696 /// Get a pointer to the derivative of the mapping from face to bulk
1697 /// coordinates.
1699 const int& face_index) const
1700 {
1701 if (face_index == 1)
1702 {
1704 }
1705 else if (face_index == -1)
1706 {
1708 }
1709 else if (face_index == -2)
1710 {
1712 }
1713 else if (face_index == 2)
1714 {
1716 }
1717 else if (face_index == -3)
1718 {
1720 }
1721 else if (face_index == 3)
1722 {
1724 }
1725 else
1726 {
1727 std::string err = "Face index should be in {-1, +1}.";
1728 throw OomphLibError(
1729 err, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
1730 }
1731 }
1732 };
1733
1734 //=======================================================================
1735 /// SolidQElement elements are quadrilateral elements whose
1736 /// derivatives also include those based upon the lagrangian
1737 /// positions of the nodes.
1738 /// They are the basis for solid mechanics elements
1739 //=======================================================================
1740 template<unsigned DIM, unsigned NNODE_1D>
1742 {
1743 };
1744
1745
1746 //=======================================================================
1747 /// SolidQElement elements, specialised to one spatial dimension
1748 //=======================================================================
1749 template<unsigned NNODE_1D>
1750 class SolidQElement<1, NNODE_1D> : public virtual QElement<1, NNODE_1D>,
1751 public virtual QSolidElementBase
1752 {
1753 public:
1754 /// Constructor
1756 {
1757 // Set the Lagrangian dimension of the element
1758 set_lagrangian_dimension(1);
1759 }
1760
1761 /// Broken copy constructor
1763
1764 /// Broken assignment operator
1765 /*void operator=(const SolidQElement&) = delete;*/
1766
1767 /// Overload the output function
1768 void output(std::ostream& outfile)
1769 {
1770 FiniteElement::output(outfile);
1771 }
1772
1773 /// Output at n_plot points
1774 inline void output(std::ostream& outfile, const unsigned& n_plot);
1775
1776 /// C-style output
1777 void output(FILE* file_pt)
1778 {
1779 FiniteElement::output(file_pt);
1780 }
1781
1782 /// C_style output at n_plot points
1783 inline void output(FILE* file_pt, const unsigned& n_plot);
1784
1785 /// Build the lower-dimensional FaceElement (an element of type
1786 /// SolidQElement<0,NNODE_1D>).
1787 /// The face index takes one of two values corresponding
1788 /// to the two possible faces:
1789 /// -1 (Left) s[0] = -1.0
1790 /// +1 (Right) s[0] = 1.0
1791 inline void build_face_element(const int& face_index,
1792 FaceElement* face_element_pt);
1793 };
1794
1795 // For the dumb Intel 9.0 compiler, these need to live in here
1796
1797 /// ////////////////////////////////////////////////////////////////////////
1798 /// ////////////////////////////////////////////////////////////////////////
1799 // SolidQElements
1800 /// ////////////////////////////////////////////////////////////////////////
1801 /// ////////////////////////////////////////////////////////////////////////
1802
1803 /// ////////////////////////////////////////////////////////////////////////
1804 // 1D SolidQElements
1805 /// ////////////////////////////////////////////////////////////////////////
1806
1807
1808 //=======================================================================
1809 /// The output function for n_plot points in each coordinate direction
1810 /// for the 1D element
1811 //=======================================================================
1812 template<unsigned NNODE_1D>
1813 void SolidQElement<1, NNODE_1D>::output(std::ostream& outfile,
1814 const unsigned& n_plot)
1815 {
1816 // Local variables
1817 Vector<double> s(1);
1818
1819 // Tecplot header info
1820 outfile << "ZONE I=" << n_plot << std::endl;
1821
1822 // Find the dimension of the nodes
1823 unsigned n_dim = this->nodal_dimension();
1824
1825 // Find the Lagrangian dimension of the first node
1826 unsigned n_lagr = static_cast<SolidNode*>(node_pt(0))->nlagrangian();
1827
1828 // Loop over element nodes
1829 for (unsigned l = 0; l < n_plot; l++)
1830 {
1831 s[0] = -1.0 + l * 2.0 / (n_plot - 1);
1832
1833 // Output the Eulerian coordinates
1834 for (unsigned i = 0; i < n_dim; i++)
1835 {
1836 outfile << QElement<1, NNODE_1D>::interpolated_x(s, i) << " ";
1837 }
1838 // Output the Lagranian coordinates
1839 for (unsigned i = 0; i < n_lagr; i++)
1840 {
1841 outfile << SolidQElement<1, NNODE_1D>::interpolated_xi(s, i) << " ";
1842 }
1843 outfile << std::endl;
1844 }
1845 outfile << std::endl;
1846 }
1847
1848 //=======================================================================
1849 /// C-style output function for n_plot points in each coordinate direction
1850 /// for the 1D element
1851 //=======================================================================
1852 template<unsigned NNODE_1D>
1853 void SolidQElement<1, NNODE_1D>::output(FILE* file_pt, const unsigned& n_plot)
1854 {
1855 // Local variables
1856 Vector<double> s(1);
1857
1858 // Tecplot header info
1859 fprintf(file_pt, "ZONE I=%i\n", n_plot);
1860
1861 // Find the dimension of the nodes
1862 unsigned n_dim = this->nodal_dimension();
1863
1864 // Find the Lagrangian dimension of the first node
1865 unsigned n_lagr = static_cast<SolidNode*>(node_pt(0))->nlagrangian();
1866
1867 // Loop over element nodes
1868 for (unsigned l = 0; l < n_plot; l++)
1869 {
1870 s[0] = -1.0 + l * 2.0 / (n_plot - 1);
1871
1872 // Output the Eulerian coordinates
1873 for (unsigned i = 0; i < n_dim; i++)
1874 {
1875 fprintf(file_pt, "%g ", QElement<1, NNODE_1D>::interpolated_x(s, i));
1876 }
1877 // Output the Lagranian coordinates
1878 for (unsigned i = 0; i < n_lagr; i++)
1879 {
1880 fprintf(
1882 }
1883 fprintf(file_pt, "\n");
1884 }
1885 fprintf(file_pt, "\n");
1886 }
1887
1888
1889 //===========================================================
1890 /// Function to setup geometrical information for lower-dimensional
1891 /// FaceElements (which are of type SolidQElement<0,1>).
1892 //===========================================================
1893 template<unsigned NNODE_1D>
1895 const int& face_index, FaceElement* face_element_pt)
1896 {
1897 // Build the standard non-solid FaceElement
1898 QElement<1, NNODE_1D>::build_face_element(face_index, face_element_pt);
1899
1900 // Set the Lagrangian dimension from the first node of the present element
1901 dynamic_cast<SolidFiniteElement*>(face_element_pt)
1902 ->set_lagrangian_dimension(
1903 static_cast<SolidNode*>(node_pt(0))->nlagrangian());
1904 }
1905
1906
1907 //=======================================================================
1908 /// SolidQElement elements, specialised to two spatial dimensions
1909 //=======================================================================
1910 template<unsigned NNODE_1D>
1911 class SolidQElement<2, NNODE_1D> : public virtual QElement<2, NNODE_1D>,
1912 public virtual QSolidElementBase
1913 {
1914 public:
1915 /// Constructor
1917 : QElementBase(),
1918 QElement<2, NNODE_1D>(),
1921 {
1922 // Set the Lagrangian dimension of the element
1923 set_lagrangian_dimension(2);
1924 }
1925
1926 /// Broken copy constructor
1928
1929 /// Broken assignment operator
1930 /*void operator=(const SolidQElement&) = delete;*/
1931
1932 /// Overload the output function
1933 void output(std::ostream& outfile)
1934 {
1935 FiniteElement::output(outfile);
1936 }
1937
1938 /// Output at n_plot^2 points
1939 inline void output(std::ostream& outfile, const unsigned& n_plot);
1940
1941 /// C-style output
1942 void output(FILE* file_pt)
1943 {
1944 FiniteElement::output(file_pt);
1945 }
1946
1947 /// C_style output at n_plot points
1948 inline void output(FILE* file_pt, const unsigned& n_plot);
1949
1950
1951 /// Build the lower-dimensional FaceElement (an element of type
1952 /// SolidQElement<1,NNODE_1D>).The face index takes one of four values
1953 /// corresponding to the four possible faces:
1954 /// -1 (West) s[0] = -1.0
1955 /// +1 (East) s[0] = 1.0
1956 /// -2 (South) s[1] = -1.0
1957 /// +2 (North) s[1] = 1.0
1958 inline void build_face_element(const int& face_index,
1959 FaceElement* face_element_pt);
1960 };
1961
1962
1963 /// ////////////////////////////////////////////////////////////////////////
1964 // 2D SolidQElements
1965 /// ////////////////////////////////////////////////////////////////////////
1966
1967 //===========================================================
1968 /// The output function for any number of points per element
1969 //===========================================================
1970 template<unsigned NNODE_1D>
1971 void SolidQElement<2, NNODE_1D>::output(std::ostream& outfile,
1972 const unsigned& n_plot)
1973 {
1974 // Local variables
1975 Vector<double> s(2);
1976
1977 // Tecplot header info
1978 outfile << "ZONE I=" << n_plot << ", J=" << n_plot << std::endl;
1979
1980 // Find the dimension of the nodes
1981 unsigned n_dim = this->nodal_dimension();
1982
1983 // Find the Lagrangian dimension of the first node
1984 unsigned n_lagr = static_cast<SolidNode*>(node_pt(0))->nlagrangian();
1985
1986 // Loop over plot points
1987 for (unsigned l2 = 0; l2 < n_plot; l2++)
1988 {
1989 s[1] = -1.0 + l2 * 2.0 / (n_plot - 1);
1990 for (unsigned l1 = 0; l1 < n_plot; l1++)
1991 {
1992 s[0] = -1.0 + l1 * 2.0 / (n_plot - 1);
1993
1994 // Output the Eulerian coordinates
1995 for (unsigned i = 0; i < n_dim; i++)
1996 {
1997 outfile << QElement<2, NNODE_1D>::interpolated_x(s, i) << " ";
1998 }
1999 // Output the Lagranian coordinates
2000 for (unsigned i = 0; i < n_lagr; i++)
2001 {
2002 outfile << SolidQElement<2, NNODE_1D>::interpolated_xi(s, i) << " ";
2003 }
2004 outfile << std::endl;
2005 }
2006 }
2007 outfile << std::endl;
2008 }
2009
2010
2011 //====================================================================
2012 /// C-style output function for any number of points per element
2013 //====================================================================
2014 template<unsigned NNODE_1D>
2015 void SolidQElement<2, NNODE_1D>::output(FILE* file_pt, const unsigned& n_plot)
2016 {
2017 // Local variables
2018 Vector<double> s(2);
2019
2020 // Tecplot header info
2021 fprintf(file_pt, "ZONE I=%i, J=%i\n", n_plot, n_plot);
2022
2023 // Find the dimension of the nodes
2024 unsigned n_dim = this->nodal_dimension();
2025
2026 // Find the Lagrangian dimension of the first node
2027 unsigned n_lagr = static_cast<SolidNode*>(node_pt(0))->nlagrangian();
2028
2029 // Loop over plot points
2030 for (unsigned l2 = 0; l2 < n_plot; l2++)
2031 {
2032 s[1] = -1.0 + l2 * 2.0 / (n_plot - 1);
2033 for (unsigned l1 = 0; l1 < n_plot; l1++)
2034 {
2035 s[0] = -1.0 + l1 * 2.0 / (n_plot - 1);
2036
2037 // Output the Eulerian coordinates
2038 for (unsigned i = 0; i < n_dim; i++)
2039 {
2040 fprintf(file_pt, "%g ", QElement<2, NNODE_1D>::interpolated_x(s, i));
2041 }
2042 // Output the Lagranian coordinates
2043 for (unsigned i = 0; i < n_lagr; i++)
2044 {
2045 fprintf(
2047 }
2048 fprintf(file_pt, "\n");
2049 }
2050 }
2051 fprintf(file_pt, "\n");
2052 }
2053
2054
2055 //===========================================================
2056 /// Function to setup geometrical information for lower-dimensional
2057 /// FaceElements (which are of type SolidQElement<1,NNODE_1D>).
2058 //===========================================================
2059 template<unsigned NNODE_1D>
2061 const int& face_index, FaceElement* face_element_pt)
2062 {
2063 // Build the standard non-solid FaceElement
2064 QElement<2, NNODE_1D>::build_face_element(face_index, face_element_pt);
2065
2066 // Set the Lagrangian dimension from the first node of the present element
2067 dynamic_cast<SolidFiniteElement*>(face_element_pt)
2068 ->set_lagrangian_dimension(
2069 static_cast<SolidNode*>(node_pt(0))->nlagrangian());
2070 }
2071
2072
2073 //=======================================================================
2074 /// SolidQElement elements, specialised to three spatial dimensions
2075 //=======================================================================
2076 template<unsigned NNODE_1D>
2077 class SolidQElement<3, NNODE_1D> : public virtual QElement<3, NNODE_1D>,
2078 public virtual QSolidElementBase
2079 {
2080 public:
2081 /// Constructor
2083 : QElementBase(),
2084 QElement<3, NNODE_1D>(),
2087 {
2088 // Set the Lagrangian dimension of the element
2089 set_lagrangian_dimension(3);
2090 }
2091
2092 /// Broken copy constructor
2094
2095 /// Broken assignment operator
2096 /*void operator=(const SolidQElement&) = delete;*/
2097
2098 /// Overload the output function
2099 void output(std::ostream& outfile)
2100 {
2101 FiniteElement::output(outfile);
2102 }
2103
2104 /// Output at n_plot^2 points
2105 inline void output(std::ostream& outfile, const unsigned& n_plot);
2106
2107 /// C-style output
2108 void output(FILE* file_pt)
2109 {
2110 FiniteElement::output(file_pt);
2111 }
2112
2113 /// C_style output at n_plot points
2114 inline void output(FILE* file_pt, const unsigned& n_plot);
2115
2116
2117 /// Build the lower-dimensional FaceElement (an element of type
2118 /// SolidQElement<2,NNODE_1D>). The face index takes of one
2119 /// six values corresponding
2120 /// to the six possible faces:
2121 /// -1 (Left) s[0] = -1.0
2122 /// +1 (Right) s[0] = 1.0
2123 /// -2 (Down) s[1] = -1.0
2124 /// +2 (Up) s[1] = 1.0
2125 /// -3 (Back) s[2] = -1.0
2126 /// +3 (Front) s[2] = 1.0
2127 inline void build_face_element(const int& face_index,
2128 FaceElement* face_element_pt);
2129 };
2130
2131
2132 /// ////////////////////////////////////////////////////////////////////////
2133 // 3D SolidQElements
2134 /// ////////////////////////////////////////////////////////////////////////
2135
2136 //===========================================================
2137 /// The output function for any number of points per element
2138 //===========================================================
2139 template<unsigned NNODE_1D>
2140 void SolidQElement<3, NNODE_1D>::output(std::ostream& outfile,
2141 const unsigned& n_plot)
2142 {
2143 // Local variables
2144 Vector<double> s(3);
2145
2146 // Tecplot header info
2147 outfile << "ZONE I=" << n_plot << ", J=" << n_plot << ", K=" << n_plot
2148 << std::endl;
2149
2150 // Find the dimension of the nodes
2151 unsigned n_dim = this->nodal_dimension();
2152
2153 // Find the Lagrangian dimension of the first node
2154 unsigned n_lagr = static_cast<SolidNode*>(node_pt(0))->nlagrangian();
2155
2156 // Loop over plot points
2157 for (unsigned l3 = 0; l3 < n_plot; l3++)
2158 {
2159 s[2] = -1.0 + l3 * 2.0 / (n_plot - 1);
2160 for (unsigned l2 = 0; l2 < n_plot; l2++)
2161 {
2162 s[1] = -1.0 + l2 * 2.0 / (n_plot - 1);
2163 for (unsigned l1 = 0; l1 < n_plot; l1++)
2164 {
2165 s[0] = -1.0 + l1 * 2.0 / (n_plot - 1);
2166
2167 // Output the Eulerian coordinates
2168 for (unsigned i = 0; i < n_dim; i++)
2169 {
2170 outfile << QElement<3, NNODE_1D>::interpolated_x(s, i) << " ";
2171 }
2172 // Output the Lagranian coordinates
2173 for (unsigned i = 0; i < n_lagr; i++)
2174 {
2175 outfile << SolidQElement<3, NNODE_1D>::interpolated_xi(s, i) << " ";
2176 }
2177 outfile << std::endl;
2178 }
2179 }
2180 }
2181 outfile << std::endl;
2182 }
2183
2184
2185 //====================================================================
2186 /// C-style output function for any number of points per element
2187 //====================================================================
2188 template<unsigned NNODE_1D>
2189 void SolidQElement<3, NNODE_1D>::output(FILE* file_pt, const unsigned& n_plot)
2190 {
2191 // Local variables
2192 Vector<double> s(3);
2193
2194 // Tecplot header info
2195 fprintf(file_pt, "ZONE I=%i, J=%i, K=%i\n", n_plot, n_plot, n_plot);
2196
2197 // Find the dimension of the nodes
2198 unsigned n_dim = this->nodal_dimension();
2199
2200 // Find the Lagrangian dimension of the first node
2201 unsigned n_lagr = static_cast<SolidNode*>(node_pt(0))->nlagrangian();
2202
2203 // Loop over plot points
2204 for (unsigned l3 = 0; l3 < n_plot; l3++)
2205 {
2206 s[2] = -1.0 + l3 * 2.0 / (n_plot - 1);
2207 for (unsigned l2 = 0; l2 < n_plot; l2++)
2208 {
2209 s[1] = -1.0 + l2 * 2.0 / (n_plot - 1);
2210 for (unsigned l1 = 0; l1 < n_plot; l1++)
2211 {
2212 s[0] = -1.0 + l1 * 2.0 / (n_plot - 1);
2213
2214 // Output the Eulerian coordinates
2215 for (unsigned i = 0; i < n_dim; i++)
2216 {
2217 fprintf(
2218 file_pt, "%g ", QElement<3, NNODE_1D>::interpolated_x(s, i));
2219 }
2220 // Output the Lagranian coordinates
2221 for (unsigned i = 0; i < n_lagr; i++)
2222 {
2223 fprintf(file_pt,
2224 "%g ",
2226 }
2227 fprintf(file_pt, "\n");
2228 }
2229 }
2230 }
2231 fprintf(file_pt, "\n");
2232 }
2233
2234
2235 //===========================================================
2236 /// Function to setup geometrical information for lower-dimensional
2237 /// FaceElements (which are of type SolidQElement<1,NNODE_1D>).
2238 //===========================================================
2239 template<unsigned NNODE_1D>
2241 const int& face_index, FaceElement* face_element_pt)
2242 {
2243 // Build the standard non-solid FaceElement
2244 QElement<3, NNODE_1D>::build_face_element(face_index, face_element_pt);
2245
2246 // Set the Lagrangian dimension from the first node of the present element
2247 dynamic_cast<SolidFiniteElement*>(face_element_pt)
2248 ->set_lagrangian_dimension(
2249 static_cast<SolidNode*>(node_pt(0))->nlagrangian());
2250 }
2251
2252
2253 //==============================================================
2254 /// A class that is used to template the refineable Q elements
2255 /// by dimension. It's really nothing more than a policy class
2256 //=============================================================
2257 template<unsigned DIM>
2259 {
2260 public:
2261 /// Empty constuctor
2263 };
2264
2265 //==============================================================
2266 /// A class that is used to template the p-refineable Q elements
2267 /// by dimension. It's really nothing more than a policy class.
2268 /// The default template parameter ensures that these elements
2269 /// inherit from the QElement of the correct type if they start
2270 /// with a p-order higher than linear (e.g. Navier-Stokes Elements).
2271 //=============================================================
2272 template<unsigned DIM, unsigned INITIAL_NNODE_1D = 2>
2274 {
2275 public:
2276 /// Empty constuctor
2278 };
2279
2280 //==============================================================
2281 /// A class that is used to template the solid refineable Q elements
2282 /// by dimension. It's really nothing more than a policy class
2283 //=============================================================
2284 template<unsigned DIM>
2286 {
2287 public:
2288 /// Empty constuctor
2290 };
2291
2292
2293} // namespace oomph
2294
2295#endif
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
char t
Definition: cfortran.h:568
Base class for all brick elements.
Definition: Qelements.h:1233
virtual Node * vertex_node_pt(const unsigned &j) const =0
Pointer to the j-th vertex node in the element.
virtual unsigned nvertex_node() const =0
Number of vertex nodes in the element.
BrickElementBase()
Constructor. Empty.
Definition: Qelements.h:1236
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
Definition: shape.h:278
FaceElements are elements that coincide with the faces of higher-dimensional "bulk" elements....
Definition: elements.h:4338
A general Finite Element class.
Definition: elements.h:1313
virtual void set_macro_elem_pt(MacroElement *macro_elem_pt)
Set pointer to macro element – can be overloaded in derived elements to perform additional tasks.
Definition: elements.h:1872
virtual void output(std::ostream &outfile)
Output the element data — typically the values at the nodes in a format suitable for post-processing.
Definition: elements.h:3050
virtual double s_min() const
Min value of local coordinate.
Definition: elements.h:2793
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
Definition: elements.cc:3962
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition: elements.h:2611
MacroElement * Macro_elem_pt
Pointer to the element's macro element (NULL by default)
Definition: elements.h:1683
virtual double s_max() const
Max. value of local coordinate.
Definition: elements.h:2803
MacroElement * macro_elem_pt()
Access function to pointer to macro element.
Definition: elements.h:1878
virtual unsigned nnode_1d() const
Return the number of nodes along one edge of the element Default is to return zero — must be overload...
Definition: elements.h:2218
Class for multidimensional Gaussian integration rules.
Definition: integral.h:145
Base class for all line elements.
Definition: Qelements.h:466
LineElementBase()
Constructor. Empty.
Definition: Qelements.h:469
virtual Node * vertex_node_pt(const unsigned &j) const =0
Pointer to the j-th vertex node in the element.
virtual unsigned nvertex_node() const =0
Number of vertex nodes in the element.
Base class for MacroElement s that are used during mesh refinement in domains with curvlinear and/or ...
Definition: macro_element.h:73
void macro_map(const Vector< double > &s, Vector< double > &r)
The mapping from local to global coordinates at the current time : r(s)
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
An OomphLibError object which should be thrown when an run-time error is encountered....
A class that is used to template the p-refineable Q elements by dimension. It's really nothing more t...
Definition: Qelements.h:2274
PRefineableQElement()
Empty constuctor.
Definition: Qelements.h:2277
/////////////////////////////////////////////////////////////////// /////////////////////////////////...
Definition: Qelements.h:91
virtual void set_macro_elem_pt(MacroElement *macro_elem_pt)
Set pointer to macro element also sets up storage for the reference coordinates and initialises them.
Definition: Qelements.h:143
double & s_macro_ur(const unsigned &i)
Access fct to the i-th coordinate of the element's "upper right" vertex in the associated MacroElemen...
Definition: Qelements.h:202
QElementBase(const QElementBase &)=delete
Broken copy constructor.
Vector< double > * S_macro_ur_pt
Pointer to vector of upper right vertex coords. in macro element.
Definition: Qelements.h:317
QElementBase()
Constructor: Initialise pointers to macro element reference coords.
Definition: Qelements.h:94
bool local_coord_is_valid(const Vector< double > &s)
Check whether the local coordinate are valid or not.
Definition: Qelements.h:113
ElementGeometry::ElementGeometry element_geometry() const
It's a Q element!
Definition: Qelements.h:307
void get_x_from_macro_element(const unsigned &t, const Vector< double > &s, Vector< double > &x)
Global coordinates as function of local coordinates at previous time "level" t (t=0: present; t>0: pr...
Definition: Qelements.h:274
double & s_macro_ll(const unsigned &i)
Access fct to the i-th coordinate of the element's "lower left" vertex in the associated MacroElement...
Definition: Qelements.h:186
void get_x_from_macro_element(const Vector< double > &s, Vector< double > &x) const
Global coordinates as function of local coordinates. using the macro element representation.
Definition: Qelements.h:249
void move_local_coord_back_into_element(Vector< double > &s) const
Adjust local coordinates so that they're located inside the element.
Definition: Qelements.h:129
double s_macro_ll(const unsigned &i) const
Access fct to the i-th coordinate of the element's "lower left" vertex in the associated MacroElement...
Definition: Qelements.h:218
unsigned nnode_on_face() const
Return number of nodes on one face of the element. Always nnode_1d^(el_dim - 1).
Definition: Qelements.h:299
Vector< double > * S_macro_ll_pt
Pointer to vector of lower left vertex coords. in macro element.
Definition: Qelements.h:314
double s_macro_ur(const unsigned &i) const
Access fct to the i-th coordinate of the element's "upper right" vertex in the associated MacroElemen...
Definition: Qelements.h:234
virtual ~QElementBase()
Broken assignment operator.
Definition: Qelements.h:103
/////////////////////////////////////////////////////////////////// /////////////////////////////////...
Definition: Qelements.h:64
QElementGeometricBase()
Empty default constructor.
Definition: Qelements.h:67
QElementGeometricBase(const QElementGeometricBase &)=delete
Broken copy constructor.
std::string tecplot_zone_string(const unsigned &nplot) const
Return string for tecplot zone header (when plotting nplot points in each "coordinate direction)
Definition: Qelements.h:717
BulkCoordinateDerivativesFctPt bulk_coordinate_derivatives_fct_pt(const int &face_index) const
Get a pointer to the derivative of the mapping from face to bulk coordinates.
Definition: Qelements.h:797
void write_paraview_output_offset_information(std::ofstream &file_out, const unsigned &nplot, unsigned &counter) const
Fill in the offset information for paraview plot. Needs to be implemented for each new geometric elem...
Definition: Qelements.h:631
QElement(const QElement &)=delete
Broken copy constructor.
void local_fraction_of_node(const unsigned &j, Vector< double > &s_fraction)
Get the local fraction of node j in the element.
Definition: Qelements.h:589
unsigned get_bulk_node_number(const int &face_index, const unsigned &i) const
Get the number of the ith node on face face_index in the bulk node vector.
Definition: Qelements.h:739
unsigned nvertex_node() const
Number of vertex nodes in the element.
Definition: Qelements.h:549
unsigned nplot_points_paraview(const unsigned &nplot) const
Return the number of actual plot points for paraview plot with parameter nplot.
Definition: Qelements.h:616
static Gauss< 1, NNODE_1D > Default_integration_scheme
Default integration rule: Gaussian integration of same 'order' as the element.
Definition: Qelements.h:490
double invert_jacobian_mapping(const DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
Overload the template-free interface for the calculation of inverse jacobian matrix....
Definition: Qelements.h:530
CoordinateMappingFctPt face_to_bulk_coordinate_fct_pt(const int &face_index) const
Get a pointer to the function mapping face coordinates to bulk coordinates.
Definition: Qelements.h:776
double s_min() const
Min. value of local coordinate.
Definition: Qelements.h:537
unsigned nnode_1d() const
Number of nodes along each element edge.
Definition: Qelements.h:609
double local_one_d_fraction_of_node(const unsigned &n1d, const unsigned &i)
This function returns the local fraction of all nodes at the n-th position in a one dimensional expan...
Definition: Qelements.h:598
unsigned nsub_elements_paraview(const unsigned &nplot) const
Return the number of local sub-elements for paraview plot with parameter nplot.
Definition: Qelements.h:623
void write_paraview_type(std::ofstream &file_out, const unsigned &nplot) const
Return the paraview element type. Needs to be implemented for each new geometric element type; see ht...
Definition: Qelements.h:650
void local_coordinate_of_node(const unsigned &j, Vector< double > &s) const
Get local coordinates of node j in the element; vector sets its own size.
Definition: Qelements.h:580
unsigned nplot_points(const unsigned &nplot) const
Return total number of plot points (when plotting nplot points in each "coordinate direction)
Definition: Qelements.h:726
void get_s_plot(const unsigned &i, const unsigned &nplot, Vector< double > &s, const bool &use_equally_spaced_interior_sample_points=false) const
Get cector of local coordinates of plot point i (when plotting nplot points in each "coordinate direc...
Definition: Qelements.h:692
int face_outer_unit_normal_sign(const int &face_index) const
Get the sign of the outer unit normal on the face given by face_index.
Definition: Qelements.h:761
double s_max() const
Max. value of local coordinate.
Definition: Qelements.h:543
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
Definition: Qelements.h:555
void write_paraview_offsets(std::ofstream &file_out, const unsigned &nplot, unsigned &offset_sum) const
Return the offsets for the paraview sub-elements. Needs to be implemented for each new geometric elem...
Definition: Qelements.h:663
CoordinateMappingFctPt face_to_bulk_coordinate_fct_pt(const int &face_index) const
Get a pointer to the function mapping face coordinates to bulk coordinates.
Definition: Qelements.h:1172
double s_min() const
Min. value of local coordinate.
Definition: Qelements.h:894
unsigned nplot_points(const unsigned &nplot) const
Return total number of plot points (when plotting nplot points in each "coordinate direction)
Definition: Qelements.h:1107
void write_paraview_type(std::ofstream &file_out, const unsigned &nplot) const
Return the paraview element type. Needs to be implemented for each new geometric element type; see ht...
Definition: Qelements.h:1025
double invert_jacobian_mapping(const DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
Overload the template-free interface for the calculation of inverse jacobian matrix....
Definition: Qelements.h:887
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
Definition: Qelements.h:913
void write_paraview_offsets(std::ofstream &file_out, const unsigned &nplot, unsigned &offset_sum) const
Return the offsets for the paraview sub-elements. Needs to be implemented for each new geometric elem...
Definition: Qelements.h:1038
unsigned nvertex_node() const
Number of vertex nodes in the element.
Definition: Qelements.h:907
BulkCoordinateDerivativesFctPt bulk_coordinate_derivatives_fct_pt(const int &face_index) const
Get a pointer to the derivative of the mapping from face to bulk coordinates.
Definition: Qelements.h:1201
void get_s_plot(const unsigned &i, const unsigned &nplot, Vector< double > &s, const bool &use_equally_spaced_interior_sample_points=false) const
Get cector of local coordinates of plot point i (when plotting nplot points in each "coordinate direc...
Definition: Qelements.h:1067
void local_fraction_of_node(const unsigned &j, Vector< double > &s_fraction)
Get the local fraction of node j in the element.
Definition: Qelements.h:957
unsigned nnode_1d() const
Number of nodes along each element edge.
Definition: Qelements.h:979
unsigned nplot_points_paraview(const unsigned &nplot) const
Return the number of actual plot points for paraview plot with parameter nplot.
Definition: Qelements.h:986
void write_paraview_output_offset_information(std::ofstream &file_out, const unsigned &nplot, unsigned &counter) const
Fill in the offset information for paraview plot. Needs to be implemented for each new geometric elem...
Definition: Qelements.h:1001
static Gauss< 2, NNODE_1D > Default_integration_scheme
Default integration rule: Gaussian integration of same 'order' as the element.
Definition: Qelements.h:845
QElement(const QElement &)=delete
Broken copy constructor.
std::string tecplot_zone_string(const unsigned &nplot) const
Return string for tecplot zone header (when plotting nplot points in each "coordinate direction)
Definition: Qelements.h:1098
void local_coordinate_of_node(const unsigned &j, Vector< double > &s) const
Get local coordinates of node j in the element; vector sets its own size.
Definition: Qelements.h:945
unsigned nsub_elements_paraview(const unsigned &nplot) const
Return the number of local sub-elements for paraview plot with parameter nplot.
Definition: Qelements.h:993
double local_one_d_fraction_of_node(const unsigned &n1d, const unsigned &i)
This function returns the local fraction of ant nodes in the n-th positoin in a one dimensional expan...
Definition: Qelements.h:968
double s_max() const
Max. value of local coordinate.
Definition: Qelements.h:900
int face_outer_unit_normal_sign(const int &face_index) const
Get the sign of the outer unit normal on the face given by face_index.
Definition: Qelements.h:1152
unsigned get_bulk_node_number(const int &face_index, const unsigned &i) const
Get the number of the ith node on face face_index (in the bulk node vector).
Definition: Qelements.h:1120
unsigned nplot_points(const unsigned &nplot) const
Return total number of plot points (when plotting nplot points in each "coordinate direction)
Definition: Qelements.h:1572
double local_one_d_fraction_of_node(const unsigned &n1d, const unsigned &i)
This function returns the local fraction of any nodes in the n-th positoin in a one dimensional expan...
Definition: Qelements.h:1400
int face_outer_unit_normal_sign(const int &face_index) const
Get the sign of the outer unit normal on the face given by face_index.
Definition: Qelements.h:1625
BulkCoordinateDerivativesFctPt bulk_coordinate_derivatives_fct_pt(const int &face_index) const
Get a pointer to the derivative of the mapping from face to bulk coordinates.
Definition: Qelements.h:1698
double s_max() const
Max. value of local coordinate.
Definition: Qelements.h:1317
unsigned nnode_1d() const
Number of nodes along each element edge.
Definition: Qelements.h:1411
void get_s_plot(const unsigned &i, const unsigned &nplot, Vector< double > &s, const bool &use_equally_spaced_interior_sample_points=false) const
Get cector of local coordinates of plot point i (when plotting nplot points in each "coordinate direc...
Definition: Qelements.h:1526
double s_min() const
Min. value of local coordinate.
Definition: Qelements.h:1311
unsigned nplot_points_paraview(const unsigned &nplot) const
Return the number of actual plot points for paraview plot with parameter nplot.
Definition: Qelements.h:1418
CoordinateMappingFctPt face_to_bulk_coordinate_fct_pt(const int &face_index) const
Get a pointer to the function mapping face coordinates to bulk coordinates.
Definition: Qelements.h:1661
void local_fraction_of_node(const unsigned &j, Vector< double > &s_fraction)
Get the local fraction of node j in the element.
Definition: Qelements.h:1387
unsigned nvertex_node() const
Number of vertex nodes in the element.
Definition: Qelements.h:1323
unsigned get_bulk_node_number(const int &face_index, const unsigned &i) const
Get the number of the ith node on face face_index in the bulk node vector.
Definition: Qelements.h:1585
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
Definition: Qelements.h:1329
std::string tecplot_zone_string(const unsigned &nplot) const
Return string for tecplot zone header (when plotting nplot points in each "coordinate direction)
Definition: Qelements.h:1562
void write_paraview_offsets(std::ofstream &file_out, const unsigned &nplot, unsigned &offset_sum) const
Return the offsets for the paraview sub-elements. Needs to be implemented for each new geometric elem...
Definition: Qelements.h:1500
double invert_jacobian_mapping(const DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
Overload the template-free interface for the calculation of the inverse jacobian mapping....
Definition: Qelements.h:1304
QElement(const QElement &)=delete
Broken copy constructor.
void write_paraview_type(std::ofstream &file_out, const unsigned &nplot) const
Return the paraview element type. Needs to be implemented for each new geometric element type; see ht...
Definition: Qelements.h:1487
static Gauss< 3, NNODE_1D > Default_integration_scheme
Default integration rule: Gaussian integration of same 'order' as the element.
Definition: Qelements.h:1257
unsigned nsub_elements_paraview(const unsigned &nplot) const
Return the number of local sub-elements for paraview plot with parameter nplot.
Definition: Qelements.h:1425
void local_coordinate_of_node(const unsigned &j, Vector< double > &s) const
Get local coordinates of node j in the element; vector sets its own size.
Definition: Qelements.h:1372
void write_paraview_output_offset_information(std::ofstream &file_out, const unsigned &nplot, unsigned &counter) const
Fill in the offset information for paraview plot. Needs to be implemented for each new geometric elem...
Definition: Qelements.h:1433
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
Definition: Qelements.h:459
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
Definition: Qelements.h:331
virtual void set_macro_elem_pt(MacroElement *macro_elem_pt)
Broken assignment operator.
Definition: Qelements.h:349
void get_x_and_xi(const Vector< double > &s, Vector< double > &x_fe, Vector< double > &x, Vector< double > &xi_fe, Vector< double > &xi) const
Eulerian and Lagrangian coordinates as function of the local coordinates: The Eulerian position is re...
Definition: Qelements.h:379
QSolidElementBase()
Constructor: Empty.
Definition: Qelements.h:334
virtual void set_macro_elem_pt(MacroElement *macro_elem_pt, MacroElement *undeformed_macro_elem_pt)
Set pointers to "current" and "undeformed" MacroElements.
Definition: Qelements.h:360
QSolidElementBase(const QSolidElementBase &)=delete
Broken copy constructor.
Base class for all quad elements.
Definition: Qelements.h:821
QuadElementBase()
Constructor. Empty.
Definition: Qelements.h:824
virtual Node * vertex_node_pt(const unsigned &j) const =0
Pointer to the j-th vertex node in the element.
virtual unsigned nvertex_node() const =0
Number of vertex nodes in the element.
A class that is used to template the refineable Q elements by dimension. It's really nothing more tha...
Definition: Qelements.h:2259
RefineableQElement()
Empty constuctor.
Definition: Qelements.h:2262
A class that is used to template the solid refineable Q elements by dimension. It's really nothing mo...
Definition: Qelements.h:2286
RefineableSolidQElement()
Empty constuctor.
Definition: Qelements.h:2289
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
void set_undeformed_macro_elem_pt(MacroElement *undeformed_macro_elem_pt)
Set pointer to "undeformed" macro element. Can be overloaded in derived classes to perform additional...
Definition: elements.h:3699
MacroElement * undeformed_macro_elem_pt()
Access function to pointer to "undeformed" macro element.
Definition: elements.h:3706
MacroElement * Undeformed_macro_elem_pt
Pointer to the element's "undeformed" macro element (NULL by default)
Definition: elements.h:4076
virtual double interpolated_xi(const Vector< double > &s, const unsigned &i) const
Return i-th FE-interpolated Lagrangian coordinate xi[i] at local coordinate s.
Definition: elements.cc:7104
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(FILE *file_pt)
C-style output.
Definition: Qelements.h:1777
void output(std::ostream &outfile)
Broken assignment operator.
Definition: Qelements.h:1768
SolidQElement(const SolidQElement &)=delete
Broken copy constructor.
void output(FILE *file_pt)
C-style output.
Definition: Qelements.h:1942
SolidQElement(const SolidQElement &)=delete
Broken copy constructor.
void output(std::ostream &outfile)
Broken assignment operator.
Definition: Qelements.h:1933
void output(std::ostream &outfile)
Broken assignment operator.
Definition: Qelements.h:2099
void output(FILE *file_pt)
C-style output.
Definition: Qelements.h:2108
SolidQElement(const SolidQElement &)=delete
Broken copy constructor.
SolidQElement elements are quadrilateral elements whose derivatives also include those based upon the...
Definition: Qelements.h:1742
void output()
Doc the command line arguments.
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
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
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)
void faces2(const Vector< double > &s, DenseMatrix< double > &dsbulk_dsface, unsigned &interior_direction)
Function for the left and right faces, along which s2 is fixed.
void faces0(const Vector< double > &s, DenseMatrix< double > &dsbulk_dsface, unsigned &interior_direction)
Function for the back and front faces, along which s0 is fixed.
void faces1(const Vector< double > &s, DenseMatrix< double > &dsbulk_dsface, unsigned &interior_direction)
Function for the up and down faces, along which s1 is fixed.
void face4(const Vector< double > &s, Vector< double > &s_bulk)
The translation scheme for the up face (s1 = 1.0)
void face1(const Vector< double > &s, Vector< double > &s_bulk)
The translation scheme for the down face (s1 = -1.0)
void face2(const Vector< double > &s, Vector< double > &s_bulk)
The translation scheme for the back face (s2 = -1.0)
void face0(const Vector< double > &s, Vector< double > &s_bulk)
The translation scheme for the left face (s0 = -1.0)
void face5(const Vector< double > &s, Vector< double > &s_bulk)
The translation scheme for the front face (s2 = 1.0)
void face3(const Vector< double > &s, Vector< double > &s_bulk)
The translation scheme for the right face (s0 = 1.0)
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
void(* BulkCoordinateDerivativesFctPt)(const Vector< double > &s, DenseMatrix< double > &ds_bulk_dsface, unsigned &interior_direction)
Typedef for the function that returns the partial derivative of the local coordinates in the bulk ele...
Definition: elements.h:1290
void(* CoordinateMappingFctPt)(const Vector< double > &s, Vector< double > &s_bulk)
Typedef for the function that translates the face coordinate to the coordinate in the bulk element.
Definition: elements.h:1282