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