Qelements.cc
Go to the documentation of this file.
1 //LIC// ====================================================================
2 //LIC// This file forms part of oomph-lib, the object-oriented,
3 //LIC// multi-physics finite-element library, available
4 //LIC// at http://www.oomph-lib.org.
5 //LIC//
6 //LIC// 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 member functions for Qelements
31 
32 //Include the appropriate headers
33 #include "Qelements.h"
34 
35 namespace oomph
36 {
37 
38 ////////////////////////////////////////////////////////////////
39 // 1D Qelements
40 ////////////////////////////////////////////////////////////////
41 
42 //=======================================================================
43 /// Assign the static Default_integration_scheme
44 //=======================================================================
45 template<unsigned NNODE_1D>
47 
48 
49 //==================================================================
50 /// Return the node at the specified local coordinate
51 //==================================================================
52 template<unsigned NNODE_1D>
55 {
56  //Load the tolerance into a local variable
58  //There is one possible index.
59  Vector<int> index(1);
60 
61  // Determine the index
62  // -------------------
63 
64  // If we are at the lower limit, the index is zero
65  if(std::fabs(s[0] + 1.0) < tol)
66  {
67  index[0] = 0;
68  }
69  // If we are at the upper limit, the index is the number of nodes minus 1
70  else if(std::fabs(s[0] - 1.0) < tol)
71  {
72  index[0] = NNODE_1D-1;
73  }
74  // Otherwise, we have to calculate the index in general
75  else
76  {
77  // For uniformly spaced nodes the node number would be
78  double float_index = 0.5*(1.0 + s[0])*(NNODE_1D-1);
79  // Convert to an integer by taking the floor (rounding down)
80  index[0] = static_cast<int>(std::floor(float_index));
81  // What is the excess. This should be safe because the
82  // we have rounded down
83  double excess = float_index - index[0];
84  // If the excess is bigger than our tolerance there is no node,
85  // return null
86  // Note that we test at both lower and upper ends.
87  if((excess > tol) && ((1.0 - excess) > tol))
88  {
89  return 0;
90  }
91  // If we are at the upper end (i.e. the system has actually rounded up)
92  // we need to add one to the index
93  if((1.0 - excess) <= tol) {index[0] += 1;}
94  }
95  // If we've got here we have a node, so let's return a pointer to it
96  return node_pt(index[0]);
97 }
98 
99 
100 //=======================================================================
101 ///Shape function for specific QElement<1,..>
102 //=======================================================================
103 template <unsigned NNODE_1D>
105  const
106 {
107  //Local storage for the shape functions
108  double Psi[NNODE_1D];
109  //Call the OneDimensional Shape functions
110  OneDimLagrange::shape<NNODE_1D>(s[0],Psi);
111 
112  //Now let's loop over the nodal points in the element
113  //and copy the values back in
114  for(unsigned i=0;i<NNODE_1D;i++) {psi[i] = Psi[i];}
115 }
116 
117 //=======================================================================
118 ///Derivatives of shape functions for specific QElement<1,..>
119 //=======================================================================
120 template <unsigned NNODE_1D>
122  DShape &dpsids) const
123 {
124  //Local storage
125  double Psi[NNODE_1D];
126  double DPsi[NNODE_1D];
127  //Call the shape functions and derivatives
128  OneDimLagrange::shape<NNODE_1D>(s[0],Psi);
129  OneDimLagrange::dshape<NNODE_1D>(s[0],DPsi);
130 
131  //Loop over shape functions in element and assign to psi
132  for(unsigned l=0;l<NNODE_1D;l++)
133  {
134  psi[l] = Psi[l];
135  dpsids(l,0) = DPsi[l];
136  }
137 }
138 
139 //=======================================================================
140 /// Second derivatives of shape functions for specific QElement<1,..>:
141 /// d2psids(i,0) = \f$ d^2 \psi_j / d s^2 \f$
142 //=======================================================================
143 template <unsigned NNODE_1D>
145  DShape &dpsids, DShape &d2psids) const
146 {
147  //Local storage for the shape functions
148  double Psi[NNODE_1D];
149  double DPsi[NNODE_1D];
150  double D2Psi[NNODE_1D];
151  //Call the shape functions and derivatives
152  OneDimLagrange::shape<NNODE_1D>(s[0],Psi);
153  OneDimLagrange::dshape<NNODE_1D>(s[0],DPsi);
154  OneDimLagrange::d2shape<NNODE_1D>(s[0],D2Psi);
155 
156  //Loop over shape functions in element and assign to psi
157  for(unsigned l=0;l<NNODE_1D;l++)
158  {
159  psi[l] = Psi[l];
160  dpsids(l,0) = DPsi[l];
161  d2psids(l,0) = D2Psi[l];
162  }
163 }
164 
165 //=======================================================================
166 /// The output function for general 1D QElements
167 //=======================================================================
168 template <unsigned NNODE_1D>
169 void QElement<1,NNODE_1D>::output(std::ostream &outfile)
170 {
171  //Tecplot header info
172  outfile << "ZONE I=" << NNODE_1D << std::endl;
173 
174  //Find the dimension of the node
175  unsigned n_dim = this->nodal_dimension();
176 
177  //Loop over element nodes
178  for(unsigned l=0;l<NNODE_1D;l++)
179  {
180  //Loop over the dimensions and output the position
181  for(unsigned i=0;i<n_dim;i++)
182  {
183  outfile << node_pt(l)->x(i) << " ";
184  }
185  //Find out how many data values at the node
186  unsigned initial_nvalue = node_pt(l)->nvalue();
187  //Lopp over the data and output whether pinned or not
188  for(unsigned i=0;i<initial_nvalue;i++)
189  {
190  outfile << node_pt(l)->is_pinned(i) << " ";
191  }
192  outfile << std::endl;
193  }
194  outfile << std::endl;
195 }
196 
197 //=======================================================================
198 /// The output function for n_plot points in each coordinate direction
199 //=======================================================================
200 template <unsigned NNODE_1D>
201 void QElement<1,NNODE_1D>::output(std::ostream &outfile,
202  const unsigned &n_plot)
203 {
204  //Local variables
205  Vector<double> s(1);
206 
207  //Tecplot header info
208  outfile << "ZONE I=" << n_plot << std::endl;
209 
210  //Find the dimension of the nodes
211  unsigned n_dim = this->nodal_dimension();
212 
213  //Loop over plot points
214  for(unsigned l=0;l<n_plot;l++)
215  {
216  s[0] = -1.0 + l*2.0/(n_plot-1);
217  //Output the coordinates
218  for (unsigned i=0;i<n_dim;i++)
219  {
220  outfile << interpolated_x(s,i) << " " ;
221  }
222  outfile << std::endl;
223  }
224  outfile << std::endl;
225 }
226 
227 
228 
229 //=======================================================================
230 /// C style output function for general 1D QElements
231 //=======================================================================
232 template <unsigned NNODE_1D>
233 void QElement<1,NNODE_1D>::output(FILE* file_pt)
234 {
235 
236  //Tecplot header info
237  fprintf(file_pt,"ZONE I=%i\n",NNODE_1D);
238 
239  //Find the dimension of the nodes
240  unsigned n_dim = this->nodal_dimension();
241 
242  //Loop over element nodes
243  for(unsigned l=0;l<NNODE_1D;l++)
244  {
245  //Loop over the dimensions and output the position
246  for(unsigned i=0;i<n_dim;i++)
247  {
248  //outfile << Node_pt[l]->x(i) << " ";
249  fprintf(file_pt,"%g ",node_pt(l)->x(i));
250  }
251  //Find out how many data values at the node
252  unsigned initial_nvalue = node_pt(l)->nvalue();
253  //Lopp over the data and output whether pinned or not
254  for(unsigned i=0;i<initial_nvalue;i++)
255  {
256  //outfile << Node_pt[l]->is_pinned(i) << " ";
257  fprintf(file_pt,"%i ",node_pt(l)->is_pinned(i));
258  }
259  //outfile << std::endl;
260  fprintf(file_pt,"\n");
261  }
262  //outfile << std::endl;
263  fprintf(file_pt,"\n");
264 }
265 
266 //=======================================================================
267 /// C style output function for n_plot points in each coordinate direction
268 //=======================================================================
269 template <unsigned NNODE_1D>
270 void QElement<1,NNODE_1D>::output(FILE* file_pt, const unsigned &n_plot)
271 {
272  //Local variables
273  Vector<double> s(1);
274 
275  //Tecplot header info
276  //outfile << "ZONE I=" << n_plot << std::endl;
277  fprintf(file_pt,"ZONE I=%i\n",n_plot);
278 
279  //Find the dimension of the first node
280  unsigned n_dim = this->nodal_dimension();
281 
282  //Loop over plot points
283  for(unsigned l=0;l<n_plot;l++)
284  {
285  s[0] = -1.0 + l*2.0/(n_plot-1);
286 
287  //Output the coordinates
288  for (unsigned i=0;i<n_dim;i++)
289  {
290  //outfile << interpolated_x(s,i) << " " ;
291  fprintf(file_pt,"%g ",interpolated_x(s,i));
292  }
293  //outfile << std::endl;
294  fprintf(file_pt,"\n");
295  }
296  //outfile << std::endl;
297  fprintf(file_pt,"\n");
298 }
299 
300 
301 ////////////////////////////////////////////////////////////////
302 // 2D Qelements
303 ////////////////////////////////////////////////////////////////
304 
305 /// Assign the spatial integration scheme
306 template<unsigned NNODE_1D>
308 
309 
310 //==================================================================
311 /// Return the node at the specified local coordinate
312 //==================================================================
313 template<unsigned NNODE_1D>
316 {
317  //Load the tolerance into a local variable
319  //There are two possible indices.
320  Vector<int> index(2);
321  //Loop over the coordinates and determine the indices
322  for(unsigned i=0;i<2;i++)
323  {
324  //If we are at the lower limit, the index is zero
325  if(std::fabs(s[i] + 1.0) < tol)
326  {
327  index[i] = 0;
328  }
329  //If we are at the upper limit, the index is the number of nodes minus 1
330  else if(std::fabs(s[i] - 1.0) < tol)
331  {
332  index[i] = NNODE_1D-1;
333  }
334  //Otherwise, we have to calculate the index in general
335  else
336  {
337  //For uniformly spaced nodes the node number would be
338  double float_index = 0.5*(1.0 + s[i])*(NNODE_1D-1);
339  //Convert to an integer by taking the floor (rounding down)
340  index[i] = static_cast<int>(std::floor(float_index));
341  //What is the excess. This should be safe because the
342  //we have rounded down
343  double excess = float_index - index[i];
344  //If the excess is bigger than our tolerance there is no node,
345  //return null
346  //Note that we test at both lower and upper ends.
347  if((excess > tol) && ((1.0 - excess) > tol))
348  {
349  return 0;
350  }
351  //If we are at the upper end (i.e. the system has actually rounded up)
352  //we need to add one to the index
353  if((1.0 - excess) <= tol) {index[i] += 1;}
354  }
355  }
356  //If we've got here we have a node, so let's return a pointer to it
357  return node_pt(index[0] + NNODE_1D*index[1]);
358 }
359 
360 
361 
362 //=======================================================================
363 /// Shape function for specific QElement<2,..>
364 ///
365 //=======================================================================
366 template <unsigned NNODE_1D>
368  const
369 {
370  //Local storage
371  double Psi[2][NNODE_1D];
372  //Call the OneDimensional Shape functions
373  OneDimLagrange::shape<NNODE_1D>(s[0],Psi[0]);
374  OneDimLagrange::shape<NNODE_1D>(s[1],Psi[1]);
375  //Index for the shape functions
376  unsigned index=0;
377 
378  //Now let's loop over the nodal points in the element
379  //s1 is the "x" coordinate, s2 the "y"
380  for(unsigned i=0;i<NNODE_1D;i++)
381  {
382  for(unsigned j=0;j<NNODE_1D;j++)
383  {
384  //Multiply the two 1D functions together to get the 2D function
385  psi[index] = Psi[1][i]*Psi[0][j];
386  //Incremenet the index
387  ++index;
388  }
389  }
390 }
391 
392 //=======================================================================
393 ///Derivatives of shape functions for specific QElement<2,..>
394 //=======================================================================
395 template <unsigned NNODE_1D>
397  DShape &dpsids) const
398 {
399  //Local storage
400  double Psi[2][NNODE_1D];
401  double DPsi[2][NNODE_1D];
402  unsigned index=0;
403 
404  //Call the shape functions and derivatives
405  OneDimLagrange::shape<NNODE_1D>(s[0],Psi[0]);
406  OneDimLagrange::shape<NNODE_1D>(s[1],Psi[1]);
407  OneDimLagrange::dshape<NNODE_1D>(s[0],DPsi[0]);
408  OneDimLagrange::dshape<NNODE_1D>(s[1],DPsi[1]);
409 
410  //Loop over shape functions in element
411  for(unsigned i=0;i<NNODE_1D;i++)
412  {
413  for(unsigned j=0;j<NNODE_1D;j++)
414  {
415  //Assign the values
416  dpsids(index,0) = Psi[1][i]*DPsi[0][j];
417  dpsids(index,1) = DPsi[1][i]* Psi[0][j];
418  psi[index] = Psi[1][i]* Psi[0][j];
419  //Increment the index
420  ++index;
421  }
422  }
423 }
424 
425 
426 //=======================================================================
427 /// Second derivatives of shape functions for specific QElement<2,..>:
428 /// d2psids(i,0) = \f$ \partial^2 \psi_j / \partial s_0^2 \f$
429 /// d2psids(i,1) = \f$ \partial^2 \psi_j / \partial s_1^2 \f$
430 /// d2psids(i,2) = \f$ \partial^2 \psi_j / \partial s_0 \partial s_1 \f$
431 //=======================================================================
432 template <unsigned NNODE_1D>
434  DShape &dpsids, DShape &d2psids) const
435 {
436  //Local storage
437  double Psi[2][NNODE_1D];
438  double DPsi[2][NNODE_1D];
439  double D2Psi[2][NNODE_1D];
440  //Index for the assembly process
441  unsigned index=0;
442 
443  //Call the shape functions and derivatives
444  OneDimLagrange::shape<NNODE_1D>(s[0],Psi[0]);
445  OneDimLagrange::shape<NNODE_1D>(s[1],Psi[1]);
446  OneDimLagrange::dshape<NNODE_1D>(s[0],DPsi[0]);
447  OneDimLagrange::dshape<NNODE_1D>(s[1],DPsi[1]);
448  OneDimLagrange::d2shape<NNODE_1D>(s[0],D2Psi[0]);
449  OneDimLagrange::d2shape<NNODE_1D>(s[1],D2Psi[1]);
450 
451  //Loop over shape functions in element
452  for(unsigned i=0;i<NNODE_1D;i++)
453  {
454  for(unsigned j=0;j<NNODE_1D;j++)
455  {
456  //Assign the values
457  psi[index] = Psi[1][i]*Psi[0][j];
458  //First derivatives
459  dpsids(index,0) = Psi[1][i]*DPsi[0][j];
460  dpsids(index,1) = DPsi[1][i]* Psi[0][j];
461  //Second derivatives
462  //N.B. index 2 is the mixed derivative
463  d2psids(index,0) = Psi[1][i]*D2Psi[0][j];
464  d2psids(index,1) = D2Psi[1][i]* Psi[0][j];
465  d2psids(index,2) = DPsi[1][i]* DPsi[0][j];
466  //Increment the index
467  ++index;
468  }
469  }
470 }
471 
472 
473 
474 
475 //===========================================================
476 /// The output function for QElement<2,NNODE_1D>
477 //===========================================================
478 template <unsigned NNODE_1D>
479 void QElement<2,NNODE_1D>::output(std::ostream &outfile)
480 {
481  //Find the dimension of the nodes
482  unsigned n_dim = this->nodal_dimension();
483  //Node number
484  unsigned l=0;
485 
486  //Tecplot header info
487  outfile << "ZONE I=" << NNODE_1D << ", J=" << NNODE_1D << std::endl;
488  //Loop over element nodes
489  for(unsigned l2=0;l2<NNODE_1D;l2++)
490  {
491  for(unsigned l1=0;l1<NNODE_1D;l1++)
492  {
493  //Loop over the dimensions and output the position
494  for(unsigned i=0;i<n_dim;i++)
495  {
496  outfile << node_pt(l)->x(i) << " ";
497  }
498  //Find out how many data values at the node
499  unsigned initial_nvalue = node_pt(l)->nvalue();
500  //Loop over the data and output whether pinned or not
501  for(unsigned i=0;i<initial_nvalue;i++)
502  {
503  outfile << node_pt(l)->is_pinned(i) << " ";
504  }
505  outfile << std::endl;
506 
507  //Increase the node number
508  ++l;
509  }
510  }
511  outfile << std::endl;
512 }
513 
514 //=======================================================================
515 /// The output function for n_plot points in each coordinate direction
516 //=======================================================================
517 template <unsigned NNODE_1D>
518 void QElement<2,NNODE_1D>::output(std::ostream &outfile,
519  const unsigned &n_plot)
520 {
521  //Local variables
522  Vector<double> s(2);
523 
524  //Tecplot header info
525  outfile << "ZONE I=" << n_plot << ", J=" << n_plot << std::endl;
526 
527  //Find the dimension of the nodes
528  unsigned n_dim = this->nodal_dimension();
529 
530  //Loop over plot points
531  for(unsigned l2=0;l2<n_plot;l2++)
532  {
533  s[1] = -1.0 + l2*2.0/(n_plot-1);
534  for(unsigned l1=0;l1<n_plot;l1++)
535  {
536  s[0] = -1.0 + l1*2.0/(n_plot-1);
537 
538  //Output the coordinates
539  for (unsigned i=0;i<n_dim;i++)
540  {
541  outfile << interpolated_x(s,i) << " " ;
542  }
543  outfile << std::endl;
544  }
545  }
546  outfile << std::endl;
547 }
548 
549 
550 
551 
552 //===========================================================
553 /// C-style output function for QElement<2,NNODE_1D>
554 //===========================================================
555 template <unsigned NNODE_1D>
556 void QElement<2,NNODE_1D>::output(FILE* file_pt)
557 {
558  //Tecplot header info
559  //outfile << "ZONE I=" << NNODE_1D << ", J=" << NNODE_1D << std::endl;
560  fprintf(file_pt,"ZONE I=%i, J=%i\n",NNODE_1D,NNODE_1D);
561 
562  //Find the dimensions of the nodes
563  unsigned n_dim = this->nodal_dimension();
564  //Node number
565  unsigned l=0;
566 
567  //Loop over element nodes
568  for(unsigned l2=0;l2<NNODE_1D;l2++)
569  {
570  for(unsigned l1=0;l1<NNODE_1D;l1++)
571  {
572  //Loop over the dimensions and output the position
573  for(unsigned i=0;i<n_dim;i++)
574  {
575  //outfile << Node_pt[l]->x(i) << " ";
576  fprintf(file_pt,"%g ",node_pt(l)->x(i));
577  }
578  //Find out how many data values at the node
579  unsigned initial_nvalue = node_pt(l)->nvalue();
580  //Loop over the data and output whether pinned or not
581  for(unsigned i=0;i<initial_nvalue;i++)
582  {
583  //outfile << Node_pt[l]->is_pinned(i) << " ";
584  fprintf(file_pt,"%i ",node_pt(l)->is_pinned(i));
585  }
586  //outfile << std::endl;
587  fprintf(file_pt,"\n");
588  //Increase the node number
589  ++l;
590  }
591  }
592  //outfile << std::endl;
593  fprintf(file_pt,"\n");
594 }
595 
596 //=======================================================================
597 /// C-style output function for n_plot points in each coordinate direction
598 //=======================================================================
599 template <unsigned NNODE_1D>
600 void QElement<2,NNODE_1D>::output(FILE* file_pt, const unsigned &n_plot)
601 {
602  //Local variables
603  Vector<double> s(2);
604 
605  //Find the dimension of the nodes
606  unsigned n_dim = this->nodal_dimension();
607 
608  //Tecplot header info
609  //outfile << "ZONE I=" << n_plot << ", J=" << n_plot << std::endl;
610  fprintf(file_pt,"ZONE I=%i, J=%i\n",n_plot,n_plot);
611 
612  //Loop over element nodes
613  for(unsigned l2=0;l2<n_plot;l2++)
614  {
615  s[1] = -1.0 + l2*2.0/(n_plot-1);
616  for(unsigned l1=0;l1<n_plot;l1++)
617  {
618  s[0] = -1.0 + l1*2.0/(n_plot-1);
619 
620  //Output the coordinates
621  for (unsigned i=0;i<n_dim;i++)
622  {
623  //outfile << interpolated_x(s,i) << " " ;
624  fprintf(file_pt,"%g ",interpolated_x(s,i));
625 
626  }
627  //outfile << std::endl;
628  fprintf(file_pt,"\n");
629  }
630  }
631  //outfile << std::endl;
632  fprintf(file_pt,"\n");
633 }
634 
635 
636 
637 ////////////////////////////////////////////////////////////////
638 // 3D Qelements
639 ////////////////////////////////////////////////////////////////
640 
641 /// Assign the spatial integration scheme
642 template<unsigned NNODE_1D>
644 
645 //==================================================================
646 /// Return the node at the specified local coordinate
647 //==================================================================
648 template<unsigned NNODE_1D>
651 {
652  //Load the tolerance into a local variable
654  //There are now three possible indices
655  Vector<int> index(3);
656  //Loop over the coordinates
657  for(unsigned i=0;i<3;i++)
658  {
659  //If we are at the lower limit, the index is zero
660  if(std::fabs(s[i] + 1.0) < tol)
661  {
662  index[i] = 0;
663  }
664  //If we are at the upper limit, the index is the number of nodes minus 1
665  else if(std::fabs(s[i] - 1.0) < tol)
666  {
667  index[i] = NNODE_1D-1;
668  }
669  //Otherwise, we have to calculate the index in general
670  else
671  {
672  //For uniformly spaced nodes the node number would be
673  double float_index = 0.5*(1.0 + s[i])*(NNODE_1D-1);
674  //Conver to an integer by taking the floor (rounding down)
675  index[i] = static_cast<int>(std::floor(float_index));
676  //What is the excess. This should be safe because
677  //we have rounded down
678  double excess = float_index - index[i];
679  //If the excess is bigger than our tolerance there is no node,
680  //return null. Note that we test at both ends
681  if((excess > tol) && ((1.0 - excess) > tol))
682  {
683  return 0;
684  }
685  //If we are at the upper end (i.e. the system has actually rounded up)
686  //we need to add one to the index
687  if((1.0 - excess) <= tol) {index[i] += 1;}
688  }
689  }
690  //If we've got here we have a node, so let's return a pointer to it
691  return node_pt(index[0] + NNODE_1D*index[1] + NNODE_1D*NNODE_1D*index[2]);
692 }
693 
694 
695 
696 //=======================================================================
697 /// Shape function for specific QElement<3,..>
698 //=======================================================================
699 template <unsigned NNODE_1D>
701  const
702 {
703  //Local storage
704  double Psi[3][NNODE_1D];
705 
706  //Call the OneDimensional Shape functions
707  OneDimLagrange::shape<NNODE_1D>(s[0],Psi[0]);
708  OneDimLagrange::shape<NNODE_1D>(s[1],Psi[1]);
709  OneDimLagrange::shape<NNODE_1D>(s[2],Psi[2]);
710 
711  //Index for the shape functions
712  unsigned index=0;
713 
714  //Now let's loop over the nodal points in the element
715  //s1 is the "x" coordinate, s2 the "y"
716  for(unsigned i=0;i<NNODE_1D;i++)
717  {
718  for(unsigned j=0;j<NNODE_1D;j++)
719  {
720  for(unsigned k=0;k<NNODE_1D;k++)
721  {
722  /*Multiply the three 1D functions together to get the 3D function*/
723  psi[index] = Psi[2][i]*Psi[1][j]*Psi[0][k];
724  //Increment the index
725  ++index;
726  }
727  }
728  }
729 }
730 
731 //=======================================================================
732 /// Derivatives of shape functions for specific QElement<3,..>
733 //=======================================================================
734 template <unsigned NNODE_1D>
736  DShape &dpsids) const
737 {
738  //Local storage
739  double Psi[3][NNODE_1D];
740  double DPsi[3][NNODE_1D];
741  //Index of the total shape function
742  unsigned index=0;
743 
744  //Call the OneDimensional Shape functions and derivatives
745  OneDimLagrange::shape<NNODE_1D>(s[0],Psi[0]);
746  OneDimLagrange::shape<NNODE_1D>(s[1],Psi[1]);
747  OneDimLagrange::shape<NNODE_1D>(s[2],Psi[2]);
748  OneDimLagrange::dshape<NNODE_1D>(s[0],DPsi[0]);
749  OneDimLagrange::dshape<NNODE_1D>(s[1],DPsi[1]);
750  OneDimLagrange::dshape<NNODE_1D>(s[2],DPsi[2]);
751 
752 
753  //Loop over shape functions in element
754  for(unsigned i=0;i<NNODE_1D;i++)
755  {
756  for(unsigned j=0;j<NNODE_1D;j++)
757  {
758  for(unsigned k=0;k<NNODE_1D;k++)
759  {
760  //Assign the values
761  dpsids(index,0) = Psi[2][i] * Psi[1][j] * DPsi[0][k];
762  dpsids(index,1) = Psi[2][i] * DPsi[1][j] * Psi[0][k];
763  dpsids(index,2) = DPsi[2][i] * Psi[1][j] * Psi[0][k];
764 
765  psi[index] = Psi[2][i]*Psi[1][j]*Psi[0][k];
766  //Increment the index
767  ++index;
768  }
769  }
770  }
771 }
772 
773 //=======================================================================
774 /// Second derivatives of shape functions for specific QElement<3,..>:
775 /// d2psids(i,0) = \f$ \partial^2 \psi_j / \partial s_0^2 \f$
776 /// d2psids(i,1) = \f$ \partial^2 \psi_j / \partial s_1^2 \f$
777 /// d2psids(i,2) = \f$ \partial^2 \psi_j / \partial s_2^2 \f$
778 /// d2psids(i,3) = \f$ \partial^2 \psi_j / \partial s_0 \partial s_1 \f$
779 /// d2psids(i,4) = \f$ \partial^2 \psi_j / \partial s_0 \partial s_2 \f$
780 /// d2psids(i,5) = \f$ \partial^2 \psi_j / \partial s_1 \partial s_2 \f$
781 //=======================================================================
782 template <unsigned NNODE_1D>
784  DShape &dpsids, DShape &d2psids) const
785 {
786  //Local storage
787  double Psi[3][NNODE_1D];
788  double DPsi[3][NNODE_1D];
789  double D2Psi[3][NNODE_1D];
790  //Index of the shape function
791  unsigned index=0;
792 
793  //Call the OneDimensional Shape functions and derivatives
794  OneDimLagrange::shape<NNODE_1D>(s[0],Psi[0]);
795  OneDimLagrange::shape<NNODE_1D>(s[1],Psi[1]);
796  OneDimLagrange::shape<NNODE_1D>(s[2],Psi[2]);
797  OneDimLagrange::dshape<NNODE_1D>(s[0],DPsi[0]);
798  OneDimLagrange::dshape<NNODE_1D>(s[1],DPsi[1]);
799  OneDimLagrange::dshape<NNODE_1D>(s[2],DPsi[2]);
800  OneDimLagrange::d2shape<NNODE_1D>(s[0],D2Psi[0]);
801  OneDimLagrange::d2shape<NNODE_1D>(s[1],D2Psi[1]);
802  OneDimLagrange::d2shape<NNODE_1D>(s[2],D2Psi[2]);
803 
804  //Loop over shape functions in element
805  for(unsigned i=0;i<NNODE_1D;i++)
806  {
807  for(unsigned j=0;j<NNODE_1D;j++)
808  {
809  for(unsigned k=0;k<NNODE_1D;k++)
810  {
811  //Assign the values
812  psi[index] = Psi[2][i]*Psi[1][j]*Psi[0][k];
813 
814  dpsids(index,0) = Psi[2][i]*Psi[1][j]*DPsi[0][k];
815  dpsids(index,1) = Psi[2][i]*DPsi[1][j]*Psi[0][k];
816  dpsids(index,2) = DPsi[2][i]* Psi[1][j]*Psi[0][k];
817 
818  //Second derivative values
819  d2psids(index,0) = Psi[2][i]*Psi[1][j]*D2Psi[0][k];
820  d2psids(index,1) = Psi[2][i]*D2Psi[1][j]*Psi[0][k];
821  d2psids(index,2) = D2Psi[2][i]* Psi[1][j]*Psi[0][k];
822  //Convention for higher indices
823  //3: mixed 12, 4: mixed 13, 5: mixed 23
824  d2psids(index,3) = Psi[2][i]*DPsi[1][j]*DPsi[0][k];
825  d2psids(index,4) = DPsi[2][i]*Psi[1][j]*DPsi[0][k];
826  d2psids(index,5) = DPsi[2][i]*DPsi[1][j]*Psi[0][k];
827  //Increment the index
828  ++index;
829  }
830  }
831  }
832 }
833 
834 //===========================================================
835 /// The output function for QElement<3,NNODE_1D>
836 //===========================================================
837 template <unsigned NNODE_1D>
838 void QElement<3,NNODE_1D>::output(std::ostream &outfile)
839 {
840  //Find the dimension of the nodes
841  unsigned n_dim = this->nodal_dimension();
842  //The node number
843  unsigned l=0;
844 
845  //Tecplot header info
846  outfile << "ZONE I=" << NNODE_1D << ", J=" << NNODE_1D
847  << ", K=" << NNODE_1D<< std::endl;
848  //Loop over element nodes
849  for(unsigned l3=0;l3<NNODE_1D;l3++)
850  {
851  for(unsigned l2=0;l2<NNODE_1D;l2++)
852  {
853  for(unsigned l1=0;l1<NNODE_1D;l1++)
854  {
855 
856  //Loop over the dimensions and output the position
857  for(unsigned i=0;i<n_dim;i++)
858  {
859  outfile << node_pt(l)->x(i) << " ";
860  }
861  //Find out how many data values at the node
862  unsigned initial_nvalue = node_pt(l)->nvalue();
863  //Loop over the data and output whether pinned or not
864  for(unsigned i=0;i<initial_nvalue;i++)
865  {
866  outfile << node_pt(l)->is_pinned(i) << " ";
867  }
868  outfile << std::endl;
869  //Increase the node number
870  ++l;
871  }
872  }
873  }
874  outfile << std::endl;
875 }
876 
877 //=======================================================================
878 /// The output function for n_plot points in each coordinate direction
879 //=======================================================================
880 template <unsigned NNODE_1D>
881 void QElement<3,NNODE_1D>::output(std::ostream &outfile,
882  const unsigned &n_plot)
883 {
884  //Local variables
885  Vector<double> s(3);
886 
887  //Tecplot header info
888  outfile << "ZONE I=" << n_plot << ", J=" << n_plot
889  << ", K=" << n_plot << std::endl;
890 
891  //Find the dimension of the nodes
892  unsigned n_dim = this->nodal_dimension();
893 
894  //Loop over element nodes
895  for(unsigned l3=0;l3<n_plot;l3++)
896  {
897  s[2] = -1.0 + l3*2.0/(n_plot-1);
898  for(unsigned l2=0;l2<n_plot;l2++)
899  {
900  s[1] = -1.0 + l2*2.0/(n_plot-1);
901  for(unsigned l1=0;l1<n_plot;l1++)
902  {
903  s[0] = -1.0 + l1*2.0/(n_plot-1);
904 
905  //Output the coordinates
906  for (unsigned i=0;i<n_dim;i++)
907  {
908  outfile << interpolated_x(s,i) << " " ;
909  }
910  outfile << std::endl;
911  }
912  }
913  }
914  outfile << std::endl;
915 }
916 
917 
918 
919 
920 //===========================================================
921 /// C-style output function for QElement<3,NNODE_1D>
922 //===========================================================
923 template <unsigned NNODE_1D>
924 void QElement<3,NNODE_1D>::output(FILE* file_pt)
925 {
926  //Find the dimension of the nodes
927  unsigned n_dim = this->nodal_dimension();
928  //The node number
929  unsigned l=0;
930 
931  //Tecplot header info
932  fprintf(file_pt,"ZONE I=%i, J=%i, K=%i\n",
933  NNODE_1D,NNODE_1D,NNODE_1D);
934 
935  //Loop over element nodes
936  for(unsigned l3=0;l3<NNODE_1D;l3++)
937  {
938  for(unsigned l2=0;l2<NNODE_1D;l2++)
939  {
940  for(unsigned l1=0;l1<NNODE_1D;l1++)
941  {
942  //Loop over the dimensions and output the position
943  for(unsigned i=0;i<n_dim;i++)
944  {
945  fprintf(file_pt,"%g ",node_pt(l)->x(i));
946  }
947 
948  //Find out how many data values at the node
949  unsigned initial_nvalue = node_pt(l)->nvalue();
950 
951  //Loop over the data and output whether pinned or not
952  for(unsigned i=0;i<initial_nvalue;i++)
953  {
954  fprintf(file_pt,"%i ",node_pt(l)->is_pinned(i));
955  }
956  fprintf(file_pt,"\n");
957  //Increase the node number
958  ++l;
959  }
960  }
961  }
962  fprintf(file_pt,"\n");
963 
964 }
965 
966 //=======================================================================
967 /// C-style output function for n_plot points in each coordinate direction
968 //=======================================================================
969 template <unsigned NNODE_1D>
970 void QElement<3,NNODE_1D>::output(FILE* file_pt, const unsigned &n_plot)
971 {
972  //Local variables
973  Vector<double> s(3);
974 
975  //Tecplot header info
976  fprintf(file_pt,"ZONE I=%i, J=%i, K=%i\n",n_plot,n_plot,n_plot);
977 
978  //Find the dimension of the nodes
979  unsigned n_dim = this->nodal_dimension();
980 
981  //Loop over element nodes
982  for(unsigned l3=0;l3<n_plot;l3++)
983  {
984  s[2] = -1.0 + l3*2.0/(n_plot-1);
985  for(unsigned l2=0;l2<n_plot;l2++)
986  {
987  s[1] = -1.0 + l2*2.0/(n_plot-1);
988  for(unsigned l1=0;l1<n_plot;l1++)
989  {
990  s[0] = -1.0 + l1*2.0/(n_plot-1);
991 
992  //Output the coordinates
993  for (unsigned i=0;i<n_dim;i++)
994  {
995  fprintf(file_pt,"%g ",interpolated_x(s,i));
996  }
997  fprintf(file_pt,"\n");
998  }
999  }
1000  }
1001  fprintf(file_pt,"\n");
1002 }
1003 
1004 
1005 
1006 //===================================================================
1007 // Build required templates
1008 //===================================================================
1009 template class QElement<1,2>;
1010 template class QElement<1,3>;
1011 template class QElement<1,4>;
1012 
1013 template class QElement<2,2>;
1014 template class QElement<2,3>;
1015 template class QElement<2,4>;
1016 
1017 template class QElement<3,2>;
1018 template class QElement<3,3>;
1019 template class QElement<3,4>;
1020 
1021 }
cstr elem_len * i
Definition: cfortran.h:607
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
Node * get_node_at_local_coordinate(const Vector< double > &s) const
Get the node at the specified local coordinate.
Definition: Qelements.cc:54
static char t char * s
Definition: cfortran.h:572
static const double Node_location_tolerance
Default value for the tolerance to be used when locating nodes via local coordinates.
Definition: elements.h:1334
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.
static Gauss< 1, NNODE_1D > Default_integration_scheme
Default integration rule: Gaussian integration of same 'order' as the element.
Definition: Qelements.h:531