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