pseudosolid_node_update_elements.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: 1216 $
7 //LIC//
8 //LIC// $LastChangedDate: 2016-08-06 15:23:30 +0100 (Sat, 06 Aug 2016) $
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 //This header file contains elements that combine two element types in
31 //a generic way.
32 
33 #ifndef OOMPH_PSEUDO_SOLID_REMESH_ELEMENTS_HEADER
34 #define OOMPH_PSEUDO_SOLID_REMESH_ELEMENTS_HEADER
35 
36 #include "elements.h"
37 
38 namespace oomph
39 {
40 
41 
42 //===========================================================================
43 /// Helper namespace for pseudo-elastic elements
44 //===========================================================================
45 namespace PseudoSolidHelper
46 {
47 
48  /// \short Static variable to hold the default value for the pseudo-solid's
49  /// inertia parameter Lambda^2.
50  extern double Zero;
51 
52 }
53 
54 
55 //==========================================================================
56 /// A templated class that permits combination two different element types,
57 /// for the solution of problems in deforming domains. The first template
58 /// paremter BASIC is the standard element and the second SOLID solves
59 /// the equations that are used to control the mesh deformation.
60 //==========================================================================
61 template<class BASIC, class SOLID>
63  : public virtual BASIC, public virtual SOLID
64 
65 {
66  /// Boolean flag to indicate shape derivative method
68 
69  public:
70 
71  /// \short Constructor, call the BASIC and SOLID elements' constructors and
72  /// set the "density" parameter for solid element to zero
73  PseudoSolidNodeUpdateElement() : BASIC(), SOLID(),
75  {
76  SOLID::lambda_sq_pt()=&PseudoSolidHelper::Zero;
77  }
78 
79  /// \short Function to describe the local dofs of the element. The ostream
80  /// specifies the output stream to which the description
81  /// is written; the string stores the currently
82  /// assembled output that is ultimately written to the
83  /// output stream by Data::describe_dofs(...); it is typically
84  /// built up incrementally as we descend through the
85  /// call hierarchy of this function when called from
86  /// Problem::describe_dofs(...)
87  void describe_local_dofs(std::ostream& out,
88  const std::string& current_string) const
89  {
90  BASIC::describe_local_dofs(out,current_string);
91  SOLID::describe_local_dofs(out,current_string);
92  }
93 
94  /// \short Compute norm of solution: use the version in the BASIC
95  /// class if there's any ambiguity
96  void compute_norm(double& el_norm)
97  {
98  BASIC::compute_norm(el_norm);
99  }
100 
101  /// \short The required number of values is the sum of the two
102  unsigned required_nvalue(const unsigned &n) const
103  {return BASIC::required_nvalue(n) + SOLID::required_nvalue(n);}
104 
105  /// \short We assume that the solid stuff is stored at the end of
106  /// the nodes, i.e. its index is the number of continuously interplated
107  /// values in the BASIC equations.
109  {
110  //At the moment, we can't handle this case in generality so throw an
111  //error if the solid pressure is stored at the nodes
112  if(SOLID::solid_p_nodal_index() >= 0)
113  {
114  throw OomphLibError(
115  "Cannot handle (non-refineable) continuous solid pressure interpolation",
116  OOMPH_CURRENT_FUNCTION,
117  OOMPH_EXCEPTION_LOCATION);
118  }
119 
120  return SOLID::solid_p_nodal_index();
121  }
122 
123  /// \short Final override for the residuals function. Contributions are
124  /// added from both underlying element types
126  {
127  //Call the basic equations first
129  //Add the solid equations contribution
131  }
132 
133  /// \short Final override for jacobian function: Contributions are
134  /// included from both the underlying element types
136  DenseMatrix<double> &jacobian)
137  {
138  //Call the basic equations first
139  BASIC::fill_in_contribution_to_jacobian(residuals,jacobian);
140  //Call the solid equations
141  SOLID::fill_in_contribution_to_jacobian(residuals,jacobian);
142 
143  // Now fill in the off-diagonal entries (the shape derivatives),
144  fill_in_shape_derivatives(jacobian);
145  }
146 
147  /// \short Final override for mass matrix function: contributions
148  /// are included from both the underlying element types
150  Vector<double> &residuals,
151  DenseMatrix<double> &jacobian, DenseMatrix<double> &mass_matrix)
152  {
153  //Call the basic equations first
155  residuals,jacobian,mass_matrix);
156  //Call the solid equations
158  residuals,jacobian,mass_matrix);
159 
160  // Now fill in the off-diagonal entries (the shape derivatives),
161  fill_in_shape_derivatives(jacobian);
162  }
163 
164 
165  /// Evaluate shape derivatives by direct finite differencing
167  {Shape_derivs_by_direct_fd = true;}
168 
169  /// Evaluate shape derivatives by chain rule
171  {Shape_derivs_by_direct_fd = false;}
172 
173 
174  /// \short Fill in the shape derivatives of the BASIC equations
175  /// w.r.t. the solid position dofs
177  {
178  //Default is to use finite differences
180  {
181  this->fill_in_shape_derivatives_by_fd(jacobian);
182  }
183  //Otherwise need to do a bit more work
184  else
185  {
186  //Calculate storage requirements
187  const unsigned n_dof = this->ndof();
188  const unsigned n_node = this->nnode();
189  const unsigned nodal_dim = this->nodal_dimension();
190 
191  //If there are no nodes or dofs return
192  if((n_dof==0) || (n_node==0)) {return;}
193 
194  //Generalised dofs have NOT been considered, shout
195  if(this->nnodal_position_type()!=1)
196  {
197  throw OomphLibError(
198  "Shape derivatives do not (yet) allow for generalised position dofs\n",
199  OOMPH_CURRENT_FUNCTION,
200  OOMPH_EXCEPTION_LOCATION);
201  }
202 
203  //Storage for derivatives of residuals w.r.t. nodal coordinates
204  RankThreeTensor<double> dresidual_dnodal_coordinates(
205  n_dof,nodal_dim,n_node,0.0);
206 
207  //Get the analytic derivatives for the BASIC equations
208  BASIC::get_dresidual_dnodal_coordinates(dresidual_dnodal_coordinates);
209 
210  //Now add the appropriate contributions to the Jacobian
211  int local_unknown=0;
212 
213  //Loop over dofs
214  //(this will include the solid dofs,
215  // but all those contributions should be zero)
216  for(unsigned l=0;l<n_dof;l++)
217  {
218  //Loop over the nodes
219  for(unsigned n=0;n<n_node;n++)
220  {
221  //Loop over the position_types (only one)
222  unsigned k=0;
223  //Loop over the coordinates
224  for(unsigned i=0;i<nodal_dim;i++)
225  {
226  //Get the equation of the local unknown
227  local_unknown = this->position_local_eqn(n,k,i);
228 
229  //If not pinned, add the contribution to the Jacobian
230  if(local_unknown >= 0)
231  {
232  jacobian(l,local_unknown) +=
233  dresidual_dnodal_coordinates(l,i,n);
234  }
235  }
236  }
237  }
238  }
239  }
240 
241 
242  /// \short Fill in the derivatives of the BASIC equations
243  /// w.r.t. the solid position dofs
245  {
246 
247  // Flag to indicate if we use first or second order FD
248  // bool use_first_order_fd=false;
249 
250  //Find the number of nodes
251  const unsigned n_node = this->nnode();
252 
253  //If there aren't any nodes, then return straight away
254  if(n_node == 0) {return;}
255 
256  //Call the update function to ensure that the element is in
257  //a consistent state before finite differencing starts
258  this->update_before_solid_position_fd();
259 
260  //Get the number of position dofs and dimensions at the node
261  const unsigned n_position_type = this->nnodal_position_type();
262  const unsigned nodal_dim = this->nodal_dimension();
263 
264  //Find the number of dofs in the element
265  const unsigned n_dof = this->ndof();
266 
267  //Create residual newres vectors
268  Vector<double> residuals(n_dof);
269  Vector<double> newres(n_dof);
270  // Vector<double> newres_minus(n_dof);
271 
272  //Calculate the residuals (for the BASIC) equations
273  //Need to do this using fill_in because get_residuals will
274  //compute all residuals for the problem, which is
275  //a little ineffecient
276  for(unsigned m=0;m<n_dof;m++) {residuals[m] = 0.0;}
278 
279  //Need to determine which degrees of freedom are solid degrees of
280  //freedom
281  //A vector of booleans that will be true if the dof is associated
282  //with the solid equations
283  std::vector<bool> dof_is_solid(n_dof,false);
284 
285  //Now set all solid positional dofs in the vector
286  for(unsigned n=0;n<n_node;n++)
287  {
288  for(unsigned k=0;k<n_position_type;k++)
289  {
290  for(unsigned i=0;i<nodal_dim;i++)
291  {
292  int local_dof = this->position_local_eqn(n,k,i);
293  if(local_dof >= 0)
294  {
295  dof_is_solid[local_dof] = true;
296  }
297  }
298  }
299  }
300 
301  //Add the solid pressures (in solid elements without
302  //solid pressure the number will be zero).
303  unsigned n_solid_pres = this->npres_solid();
304  for(unsigned l=0;l<n_solid_pres;l++)
305  {
306  int local_dof = this->solid_p_local_eqn(l);
307  if(local_dof >= 0)
308  {
309  dof_is_solid[local_dof] = true;
310  }
311  }
312 
313 
314  //Integer storage for local unknown
315  int local_unknown=0;
316 
317  //Use default value defined in GeneralisedElement
318  const double fd_step = this->Default_fd_jacobian_step;
319 
320  //Loop over the nodes
321  for(unsigned n=0;n<n_node;n++)
322  {
323  //Loop over position dofs
324  for(unsigned k=0;k<n_position_type;k++)
325  {
326  //Loop over dimension
327  for(unsigned i=0;i<nodal_dim;i++)
328  {
329  //If the variable is free
330  local_unknown = this->position_local_eqn(n,k,i);
331  if(local_unknown >= 0)
332  {
333  //Store a pointer to the (generalised) Eulerian nodal position
334  double* const value_pt = &(this->node_pt(n)->x_gen(k,i));
335 
336  //Save the old value of the (generalised) Eulerian nodal position
337  const double old_var = *value_pt;
338 
339  //Increment the (generalised) Eulerian nodal position
340  *value_pt += fd_step;
341 
342  // Perform any auxialiary node updates
343  this->node_pt(n)->perform_auxiliary_node_update_fct();
344 
345  //Calculate the new residuals
346  //Need to do this using fill_in because get_residuals will
347  //compute all residuals for the problem, which is
348  //a little ineffecient
349  for(unsigned m=0;m<n_dof;m++) {newres[m] = 0.0;}
351 
352 // if (use_first_order_fd)
353  {
354  //Do forward finite differences
355  for(unsigned m=0;m<n_dof;m++)
356  {
357  //Stick the entry into the Jacobian matrix
358  //but only if it's not a solid dof
359  if(dof_is_solid[m] == false)
360  {
361  jacobian(m,local_unknown) = (newres[m] - residuals[m])/fd_step;
362  }
363  }
364  }
365 // else
366 // {
367 // //Take backwards step for the (generalised) Eulerian nodal
368 // // position
369 // node_pt(n)->x_gen(k,i) = old_var-fd_step;
370 
371 // //Calculate the new residuals at backward position
372 // //BASIC::get_residuals(newres_minus);
373 
374 // //Do central finite differences
375 // for(unsigned m=0;m<n_dof;m++)
376 // {
377 // //Stick the entry into the Jacobian matrix
378 // jacobian(m,local_unknown) =
379 // (newres[m] - newres_minus[m])/(2.0*fd_step);
380 // }
381 // }
382 
383  //Reset the (generalised) Eulerian nodal position
384  *value_pt = old_var;
385 
386  // Perform any auxialiary node updates
387  this->node_pt(n)->perform_auxiliary_node_update_fct();
388  }
389  }
390  }
391  }
392 
393  //End of finite difference loop
394  //Final reset of any slaved data
395  this->reset_after_solid_position_fd();
396  }
397 
398 
399  /// \short Specify Data that affects the geometry of the element
400  /// by adding the position Data to the set that's passed in.
401  /// (This functionality is required in FSI problems; set is used to
402  /// avoid double counting).
403  void identify_geometric_data(std::set<Data*> &geometric_data_pt)
404  {
405  //Loop over the node update data and add to the set
406  const unsigned n_node=this->nnode();
407  for(unsigned j=0;j<n_node;j++)
408  {
409  geometric_data_pt.insert(dynamic_cast<SolidNode*>(this->node_pt(j))
410  ->variable_position_pt());
411  }
412  }
413 
414 
415  ///Overload the output function: Call that of the basic element
416  void output(std::ostream &outfile) {BASIC::output(outfile);}
417 
418  /// \short Output function: Plot at n_p plot points using the basic element's
419  /// output function
420  void output(std::ostream &outfile, const unsigned &n_p)
421  {BASIC::output(outfile,n_p);}
422 
423  ///Overload the output function: Call that of the basic element
424  void output(FILE* file_pt) {BASIC::output(file_pt);}
425 
426  ///Output function is just the same as the basic equations
427  void output(FILE* file_pt, const unsigned &n_p)
428  {BASIC::output(file_pt,n_p);}
429 
430  /// \short Number of 'flux' terms for Z2 error estimation: Error estimation
431  /// is based on error in BASIC element
432  unsigned num_Z2_flux_terms() {return BASIC::num_Z2_flux_terms();}
433 
434 
435  /// \short Plot the error when compared against a given exact flux.
436  /// Also calculates the norm of the error and that of the exact flux.
437  /// Use version in BASIC element
439  std::ostream &outfile,
441  double& error, double& norm)
442  {
443  BASIC::compute_exact_Z2_error(outfile, exact_flux_pt,
444  error, norm);
445  }
446 
447  /// 'Flux' vector for Z2 error estimation: Error estimation
448  /// is based on error in BASIC element
450  {
451  BASIC::get_Z2_flux(s,flux);
452  }
453 
454  /// \short Number of vertex nodes in the element
455  unsigned nvertex_node() const
456  {return BASIC::nvertex_node();}
457 
458  /// \short Pointer to the j-th vertex node in the element
459  Node* vertex_node_pt(const unsigned& j) const
460  {return BASIC::vertex_node_pt(j);}
461 
462  /// \short Order of recovery shape functions for Z2 error estimation: Done
463  /// for BASIC element since it determines the refinement
464  unsigned nrecovery_order() {return BASIC::nrecovery_order();}
465 
466 
467  /// \short The number of "DOF types" that degrees of freedom in this element
468  /// are sub-divided into.
469  unsigned ndof_types() const
470  {
471  return BASIC::ndof_types() + SOLID::ndof_types();
472  }
473 
474  /// \short return the number of DOF types associated with the BASIC
475  /// elements in this combined element
476  unsigned nbasic_dof_types() const
477  {
478  return BASIC::ndof_types();
479  }
480 
481  /// \short return the number of DOF types associated with the SOLID
482  /// elements in this combined element
483  unsigned nsolid_dof_types() const
484  {
485  return SOLID::ndof_types();
486  }
487 
488  /// \short Create a list of pairs for all unknowns in this element,
489  /// so that the first entry in each pair contains the global equation
490  /// number of the unknown, while the second one contains the number
491  /// of the "DOF type" that this unknown is associated with.
492  /// This method combines the get_dof_numbers_for_unknowns(...)
493  /// method for the BASIC and SOLID elements. The basic elements
494  /// retain their DOF type numbering and the SOLID elements
495  /// DOF type numbers are incremented by nbasic_dof_types().
497  std::list<std::pair<unsigned long,unsigned> >& dof_lookup_list) const
498  {
499  // get the solid list
500  std::list<std::pair<unsigned long,unsigned> > solid_list;
501  SOLID::get_dof_numbers_for_unknowns(solid_list);
502 
503  // get the basic list
504  BASIC::get_dof_numbers_for_unknowns(dof_lookup_list);
505 
506  // get the number of basic dof types
507  unsigned nbasic_dof_types = BASIC::ndof_types();
508 
509  // add the solid lookup list to the basic lookup list
510  // incrementing the solid dof numbers by nbasic_dof_types
511  typedef std::list<std::pair<unsigned long,unsigned> >::iterator IT;
512  for (IT it=solid_list.begin();
513  it!=solid_list.end();it++)
514  {
515  std::pair<unsigned long,unsigned> new_pair;
516  new_pair.first = it->first;
517  new_pair.second = it->second + nbasic_dof_types;
518  dof_lookup_list.push_front(new_pair);
519  }
520  }
521 
522 };
523 
524 ///Explicit definition of the face geometry of these elements
525 template<class BASIC, class SOLID>
527 public virtual FaceGeometry<SOLID>
528 {
529  public:
530 
531  /// \short Constuctor calls the constructor of the SolidQElement
532  /// (Only the Intel compiler seems to need this!)
533  FaceGeometry() : FaceGeometry<SOLID>() {}
534 
535 };
536 
537 ///Explicit definition of the face geometry of these elements
538 template<class BASIC, class SOLID>
540 public virtual FaceGeometry<FaceGeometry<SOLID> >
541 {
542  public:
543 
544  /// \short Constuctor calls the constructor of the SolidQElement
545  /// (Only the Intel compiler seems to need this!)
547 
548  protected:
549 
550 };
551 
552 
553 
554 //===================================================================
555 /// Refineable version of the PseudoSolidNodeUpdateELement
556 //===================================================================
557 template<class BASIC, class SOLID>
558 class RefineablePseudoSolidNodeUpdateElement : public virtual BASIC,
559  public virtual SOLID
560 {
561 
562  public:
563 
564  /// \short Constructor, call the BASIC and SOLID elements' constructors and
565  /// set the "density" parameter for solid element to zero
568  BASIC(), SOLID()
569  {
570  SOLID::lambda_sq_pt()=&PseudoSolidHelper::Zero;
571  }
572 
573  /// \short Function to describe the local dofs of the element. The ostream
574  /// specifies the output stream to which the description
575  /// is written; the string stores the currently
576  /// assembled output that is ultimately written to the
577  /// output stream by Data::describe_dofs(...); it is typically
578  /// built up incrementally as we descend through the
579  /// call hierarchy of this function when called from
580  /// Problem::describe_dofs(...)
581  void describe_local_dofs(std::ostream& out,
582  const std::string& current_string) const
583  {
584  BASIC::describe_local_dofs(out,current_string);
585  SOLID::describe_local_dofs(out,current_string);
586  }
587 
588  /// \short The required number of values is the sum of the two
589  unsigned required_nvalue(const unsigned &n) const
590  {return BASIC::required_nvalue(n) + SOLID::required_nvalue(n);}
591 
592  /// \short The number of continuously interpolated values is the
593  /// sum of the SOLID and BASIC values
594  unsigned ncont_interpolated_values() const
595  {return BASIC::ncont_interpolated_values() +
596  SOLID::ncont_interpolated_values();}
597 
598  /// \short We assume that the solid stuff is stored at the end of
599  /// the nodes, i.e. its index is the number of continuously interplated
600  /// values in the BASIC equations.
602  {
603  //Find the index in the solid
604  int solid_p_index = SOLID::solid_p_nodal_index();
605  //If there is a solid pressure at the nodes, return the
606  //index after all the BASIC stuff
607  if(solid_p_index >= 0)
608  {return BASIC::ncont_interpolated_values() +
609  SOLID::solid_p_nodal_index();}
610  else {return solid_p_index;}
611  }
612 
613  ///\short Final override for residuals function: adds contributions
614  ///from both underlying element types
616  {
617  //Call the basic equations first
619  //Call the solid equations
621  }
622 
623  ///\short Final override for jacobian function: Calls get_jacobian() for
624  /// both of the underlying element types
626  DenseMatrix<double> &jacobian)
627  {
628  //Call the basic equations first
629  BASIC::fill_in_contribution_to_jacobian(residuals,jacobian);
630 
631  //Call the solid equations
632  SOLID::fill_in_contribution_to_jacobian(residuals,jacobian);
633 
634  // Now fill in the off-diagonal entries (the shape derivatives),
636  }
637 
638  /// \short Final override for mass matrix function: contributions
639  /// are included from both the underlying element types
641  Vector<double> &residuals,
642  DenseMatrix<double> &jacobian, DenseMatrix<double> &mass_matrix)
643  {
644  //Call the basic equations first
646  residuals,jacobian,mass_matrix);
647  //Call the solid equations
649  residuals,jacobian,mass_matrix);
650 
651  // Now fill in the off-diagonal entries (the shape derivatives),
653  }
654 
655 
656 
657  /// \short Fill in the derivatives of the BASIC equations
658  /// w.r.t. to the solid position dofs, taking hanging nodes
659  /// into account
661  {
662  //Find the number of nodes
663  const unsigned n_node = this->nnode();
664 
665  //If there are no nodes, return straight away
666  if(n_node == 0) {return;}
667 
668  //Call the update function to ensure that the element is in
669  //a consistent state before finite differencing starts
670  this->update_before_solid_position_fd();
671 
672 // bool use_first_order_fd=false;
673 
674  //Find the number of positional dofs and nodal dimension
675  const unsigned n_position_type = this->nnodal_position_type();
676  const unsigned nodal_dim = this->nodal_dimension();
677 
678  //Find the number of dofs in the element
679  const unsigned n_dof = this->ndof();
680 
681  //Create residual newres vectors
682  Vector<double> residuals(n_dof);
683  Vector<double> newres(n_dof);
684  // Vector<double> newres_minus(n_dof);
685 
686  //Calculate the residuals (for the BASIC) equations
687  //Need to do this using fill_in because get_residuals will
688  //compute all residuals for the problem, which is
689  //a little ineffecient
690  for(unsigned m=0;m<n_dof;m++) {residuals[m] = 0.0;}
692 
693  //Need to determine which degrees of freedom are solid degrees of
694  //freedom
695  //A vector of booleans that will be true if the dof is associated
696  //with the solid equations
697  std::vector<bool> dof_is_solid(n_dof,false);
698 
699  //Now set all solid positional dofs in the vector
700  //This is a bit more involved because we need to take account of
701  //any hanging nodes
702  for(unsigned n=0;n<n_node;n++)
703  {
704  //Get pointer to the local node
705  Node* const local_node_pt = this->node_pt(n);
706 
707  //If the node is not a hanging node
708  if(local_node_pt->is_hanging()==false)
709  {
710  for(unsigned k=0;k<n_position_type;k++)
711  {
712  for(unsigned i=0;i<nodal_dim;i++)
713  {
714  int local_dof = this->position_local_eqn(n,k,i);
715  if(local_dof >= 0)
716  {
717  dof_is_solid[local_dof] = true;
718  }
719  }
720  }
721  }
722  //Otherwise the node is hanging
723  else
724  {
725  //Find the local hanging object
726  HangInfo* hang_info_pt = local_node_pt->hanging_pt();
727  //Loop over the master nodes
728  const unsigned n_master = hang_info_pt->nmaster();
729  for(unsigned m=0;m<n_master;m++)
730  {
731  //Get the local equation numbers for the master node
732  DenseMatrix<int> Position_local_eqn_at_node
733  = this->local_position_hang_eqn(hang_info_pt->master_node_pt(m));
734 
735  //Loop over position dofs
736  for(unsigned k=0;k<n_position_type;k++)
737  {
738  //Loop over dimension
739  for(unsigned i=0;i<nodal_dim;i++)
740  {
741  int local_dof = Position_local_eqn_at_node(k,i);
742  if(local_dof >= 0)
743  {
744  dof_is_solid[local_dof] = true;
745  }
746  }
747  }
748  }
749  }
750  } //End of loop over nodes
751 
752  //Add the solid pressures (in solid elements without
753  //solid pressure the number will be zero).
754  unsigned n_solid_pres = this->npres_solid();
755  //Now is the solid pressure hanging
756  const int solid_p_index = this->solid_p_nodal_index();
757  //Find out whether the solid node is hanging
758  std::vector<bool> solid_p_is_hanging(n_solid_pres);
759  //If we have nodal solid pressures then read out the hanging status
760  if(solid_p_index >= 0)
761  {
762  //Loop over the solid dofs
763  for(unsigned l=0;l<n_solid_pres;l++)
764  {
765  solid_p_is_hanging[l] =
766  this->solid_pressure_node_pt(l)->is_hanging(solid_p_index);
767  }
768  }
769  //Otherwise the pressure is not nodal, so cannot hang
770  else
771  {
772  for(unsigned l=0;l<n_solid_pres;l++)
773  {
774  solid_p_is_hanging[l] = false;
775  }
776  }
777 
778  //Now we can loop of the dofs again to actually set that the appropriate
779  //dofs are solid
780  for(unsigned l=0;l<n_solid_pres;l++)
781  {
782  //If the solid pressure is not hanging
783  //we just read out the local equation numbers directly
784  if(solid_p_is_hanging[l] == false)
785  {
786  int local_dof = this->solid_p_local_eqn(l);
787  if(local_dof >= 0)
788  {
789  dof_is_solid[local_dof] = true;
790  }
791  }
792  //Otherwise solid pressure is hanging and we need to take
793  //care of the master nodes
794  else
795  {
796  //Find the local hanging object
797  HangInfo* hang_info_pt =
798  this->solid_pressure_node_pt(l)->hanging_pt(solid_p_index);
799  //Loop over the master nodes
800  const unsigned n_master = hang_info_pt->nmaster();
801  for(unsigned m=0;m<n_master;m++)
802  {
803  //Get the local dof
804  int local_dof = this->local_hang_eqn(
805  hang_info_pt->master_node_pt(m),solid_p_index);
806 
807  if(local_dof >= 0)
808  {
809  dof_is_solid[local_dof] = true;
810  }
811  }
812  }
813  } //end of loop over solid pressure dofs
814 
815 
816  //Used default value defined in GeneralisedElement
817  const double fd_step = this->Default_fd_jacobian_step;
818 
819  //Integer storage for local unknowns
820  int local_unknown=0;
821 
822  //Loop over the nodes
823  for(unsigned l=0;l<n_node;l++)
824  {
825  //Get the pointer to the node
826  Node* const local_node_pt = this->node_pt(l);
827 
828  //If the node is not a hanging node
829  if(local_node_pt->is_hanging()==false)
830  {
831  //Loop over position dofs
832  for(unsigned k=0;k<n_position_type;k++)
833  {
834  //Loop over dimension
835  for(unsigned i=0;i<nodal_dim;i++)
836  {
837  local_unknown = this->position_local_eqn(l,k,i);
838  //If the variable is free
839  if(local_unknown >= 0)
840  {
841  //Store a pointer to the (generalised) Eulerian nodal position
842  double* const value_pt = &(local_node_pt->x_gen(k,i));
843 
844  //Save the old value of the (generalised) Eulerian nodal position
845  const double old_var = *value_pt;
846 
847  //Increment the (generalised) Eulerian nodal position
848  *value_pt += fd_step;
849 
850  // Perform any auxialiary node updates
851  local_node_pt->perform_auxiliary_node_update_fct();
852 
853  //Calculate the new residuals
854  //Need to do this using fill_in because get_residuals will
855  //compute all residuals for the problem, which is
856  //a little ineffecient
857  for(unsigned m=0;m<n_dof;m++) {newres[m] = 0.0;}
859 
860 
861 // if (use_first_order_fd)
862  {
863  //Do forward finite differences
864  for(unsigned m=0;m<n_dof;m++)
865  {
866  //Stick the entry into the Jacobian matrix
867  //But only if it's not a solid dof
868  if(dof_is_solid[m]==false)
869  {
870  jacobian(m,local_unknown) =
871  (newres[m] - residuals[m])/fd_step;
872  }
873  }
874  }
875 // else
876 // {
877 // //Take backwards step for the (generalised) Eulerian nodal
878 // // position
879 // node_pt(l)->x_gen(k,i) = old_var-fd_step;
880 
881 // //Calculate the new residuals at backward position
882 // BASIC::get_residuals(newres_minus);
883 
884 // //Do central finite differences
885 // for(unsigned m=0;m<n_dof;m++)
886 // {
887 // //Stick the entry into the Jacobian matrix
888 // jacobian(m,local_unknown) =
889 // (newres[m] - newres_minus[m])/(2.0*fd_step);
890 // }
891 // }
892 
893  //Reset the (generalised) Eulerian nodal position
894  *value_pt = old_var;
895 
896  // Perform any auxialiary node updates
897  local_node_pt->perform_auxiliary_node_update_fct();
898 
899  }
900  }
901  }
902  }
903  //Otherwise it's a hanging node
904  else
905  {
906  //Find the local hanging object
907  HangInfo* hang_info_pt = local_node_pt->hanging_pt();
908  //Loop over the master nodes
909  const unsigned n_master = hang_info_pt->nmaster();
910  for(unsigned m=0;m<n_master;m++)
911  {
912  //Get the pointer to the master node
913  Node* const master_node_pt = hang_info_pt->master_node_pt(m);
914 
915  //Get the local equation numbers for the master node
916  DenseMatrix<int> Position_local_eqn_at_node
917  = this->local_position_hang_eqn(master_node_pt);
918 
919  //Loop over position dofs
920  for(unsigned k=0;k<n_position_type;k++)
921  {
922  //Loop over dimension
923  for(unsigned i=0;i<nodal_dim;i++)
924  {
925  local_unknown = Position_local_eqn_at_node(k,i);
926  //If the variable is free
927  if(local_unknown >= 0)
928  {
929  //Store a pointer to the (generalised) Eulerian nodal position
930  double* const value_pt = &(master_node_pt->x_gen(k,i));
931 
932  //Save the old value of the (generalised) Eulerian nodal
933  //position
934  const double old_var = *value_pt;
935 
936  //Increment the (generalised) Eulerian nodal position
937  *value_pt += fd_step;
938 
939  // Perform any auxialiary node updates
940  master_node_pt->perform_auxiliary_node_update_fct();
941 
942  //Calculate the new residuals
943  //Need to do this using fill_in because get_residuals will
944  //compute all residuals for the problem, which is
945  //a little ineffecient
946  for(unsigned m=0;m<n_dof;m++) {newres[m] = 0.0;}
948 
949 // if (use_first_order_fd)
950  {
951  //Do forward finite differences
952  for(unsigned m=0;m<n_dof;m++)
953  {
954  //Stick the entry into the Jacobian matrix
955  //But only if it's not a solid dof
956  if(dof_is_solid[m] == false)
957  {
958  jacobian(m,local_unknown) =
959  (newres[m] - residuals[m])/fd_step;
960  }
961  }
962  }
963 // else
964 // {
965 // //Take backwards step for the (generalised) Eulerian nodal
966 // // position
967 // master_node_pt->x_gen(k,i) = old_var-fd_step;
968 
969 // //Calculate the new residuals at backward position
970 // BASIC::get_residuals(newres_minus);
971 
972 // //Do central finite differences
973 // for(unsigned m=0;m<n_dof;m++)
974 // {
975 // //Stick the entry into the Jacobian matrix
976 // jacobian(m,local_unknown) =
977 // (newres[m] - newres_minus[m])/(2.0*fd_step);
978 // }
979 // }
980 
981  //Reset the (generalised) Eulerian nodal position
982  *value_pt = old_var;
983 
984  // Perform any auxialiary node updates
985  master_node_pt->perform_auxiliary_node_update_fct();
986  }
987  }
988  }
989  }
990  } //End of hanging node case
991 
992  } //End of loop over nodes
993 
994  //End of finite difference loop
995 
996  //Final reset of any slaved data
997  this->reset_after_solid_position_fd();
998  }
999 
1000 
1001  /// \short Specify Data that affects the geometry of the element
1002  /// by adding the position Data to the set that's passed in.
1003  /// (This functionality is required in FSI problems; set is used to
1004  /// avoid double counting). Refineable version includes hanging nodes
1005  void identify_geometric_data(std::set<Data*> &geometric_data_pt)
1006  {
1007  //Loop over the node update data and add to the set
1008  const unsigned n_node=this->nnode();
1009  for(unsigned j=0;j<n_node;j++)
1010  {
1011 
1012  //If the node is a hanging node
1013  if(this->node_pt(j)->is_hanging())
1014  {
1015  //Find the local hang info object
1016  HangInfo* hang_info_pt = this->node_pt(j)->hanging_pt();
1017 
1018  //Find the number of master nodes
1019  unsigned n_master = hang_info_pt->nmaster();
1020 
1021  //Loop over the master nodes
1022  for(unsigned m=0;m<n_master;m++)
1023  {
1024  //Get the m-th master node
1025  Node* Master_node_pt = hang_info_pt->master_node_pt(m);
1026 
1027  // Add to set
1028  geometric_data_pt.insert(
1029  dynamic_cast<SolidNode*>(Master_node_pt)->variable_position_pt());
1030  }
1031  }
1032  // Not hanging
1033  else
1034  {
1035  // Add node itself to set
1036  geometric_data_pt.insert(
1037  dynamic_cast<SolidNode*>(this->node_pt(j))->variable_position_pt());
1038  }
1039  }
1040  }
1041 
1042 
1043 
1044  /// \short Final override for the assign__additional_local_eqn_numbers():
1045  /// Call the version for the BASIC element
1047  {
1048  BASIC::assign_additional_local_eqn_numbers();
1049  SOLID::assign_additional_local_eqn_numbers();
1050  }
1051 
1052  /// Call rebuild_from_sons() for both of the underlying element types
1053  void rebuild_from_sons(Mesh* &mesh_pt)
1054  {
1055  BASIC::rebuild_from_sons(mesh_pt);
1056  SOLID::rebuild_from_sons(mesh_pt);
1057  }
1058 
1059  /// Call get_interpolated_values(...) for both of the underlying element types
1060  void get_interpolated_values(const unsigned& t, const Vector<double>&s,
1061  Vector<double>& values)
1062  {
1063  Vector<double> basic_values;
1064  BASIC::get_interpolated_values(t,s,basic_values);
1065  Vector<double> solid_values;
1066  SOLID::get_interpolated_values(t,s,solid_values);
1067 
1068  //Now add the basic value first
1069  for(Vector<double>::iterator it=basic_values.begin();
1070  it!=basic_values.end();++it)
1071  {
1072  values.push_back(*it);
1073  }
1074  //Then the solid
1075  for(Vector<double>::iterator it=solid_values.begin();
1076  it!=solid_values.end();++it)
1077  {
1078  values.push_back(*it);
1079  }
1080  }
1081 
1082 
1083  /// Call get_interpolated_values(...) for both of the underlying element types
1085  {
1086  Vector<double> basic_values;
1087  BASIC::get_interpolated_values(s,basic_values);
1088  Vector<double> solid_values;
1089  SOLID::get_interpolated_values(s,solid_values);
1090 
1091  //Now add the basic value first
1092  for(Vector<double>::iterator it=basic_values.begin();
1093  it!=basic_values.end();++it)
1094  {
1095  values.push_back(*it);
1096  }
1097  //Then the solid
1098  for(Vector<double>::iterator it=solid_values.begin();
1099  it!=solid_values.end();++it)
1100  {
1101  values.push_back(*it);
1102  }
1103  }
1104 
1105  /// \short We must compose the underlying interpolating nodes from
1106  /// the BASIC and SOLID equations, the BASIC ones are first
1107  Node* interpolating_node_pt(const unsigned &n,
1108  const int &value_id)
1109  {
1110  //Find the number of interpolated values in the BASIC equations
1111  int n_basic_values = BASIC::ncont_interpolated_values();
1112  //If the id is below this number, we call the BASIC functon
1113  if(value_id < n_basic_values)
1114  {return BASIC::interpolating_node_pt(n,value_id);}
1115  //Otherwise it's the solid and its value_id is the the current
1116  //it minus n_basic_values
1117  else
1118  {return SOLID::interpolating_node_pt(n,(value_id-n_basic_values));}
1119  }
1120 
1121  /// \short The pressure nodes are the corner nodes, so when value_id==0,
1122  /// the fraction is the same as the 1d node number, 0 or 1.
1124  const unsigned &i,
1125  const int &value_id)
1126  {
1127  //Find the number of interpolated values in the BASIC equations
1128  int n_basic_values = BASIC::ncont_interpolated_values();
1129  //If the id is below this number, we call the BASIC functon
1130  if(value_id < n_basic_values)
1131  {
1132  return BASIC::local_one_d_fraction_of_interpolating_node(n1d,i,value_id);
1133  }
1134  //Otherwise it's the solid and its value_id is the the current
1135  //it minus n_basic_values
1136  else
1137  {
1138  return
1139  SOLID::local_one_d_fraction_of_interpolating_node(
1140  n1d,i,(value_id-n_basic_values));
1141  }
1142  }
1143 
1144 
1145  /// \short The velocity nodes are the same as the geometric nodes. The
1146  /// pressure nodes must be calculated by using the same methods as
1147  /// the geometric nodes, but by recalling that there are only two pressure
1148  /// nodes per edge.
1150  const int &value_id)
1151  {
1152  //Find the number of interpolated values in the BASIC equations
1153  int n_basic_values = BASIC::ncont_interpolated_values();
1154  //If the id is below this number, we call the BASIC functon
1155  if(value_id < n_basic_values)
1156  {
1157  return BASIC::get_interpolating_node_at_local_coordinate(s,value_id);
1158  }
1159  //Otherwise it's the solid and its value_id is the the current
1160  //it minus n_basic_values
1161  else
1162  {
1163  return
1164  SOLID::get_interpolating_node_at_local_coordinate(
1165  s,(value_id-n_basic_values));
1166  }
1167  }
1168 
1169  /// \short The number of 1d pressure nodes is 2, otherwise we have
1170  /// the positional nodes
1171  unsigned ninterpolating_node_1d(const int &value_id)
1172  {
1173  //Find the number of interpolated values in the BASIC equations
1174  int n_basic_values = BASIC::ncont_interpolated_values();
1175  //If the id is below this number, we call the BASIC functon
1176  if(value_id < n_basic_values)
1177  {
1178  return BASIC::ninterpolating_node_1d(value_id);
1179  }
1180  //Otherwise it's the solid and its value_id is the the current
1181  //it minus n_basic_values
1182  else
1183  {
1184  return SOLID::ninterpolating_node_1d((value_id-n_basic_values));
1185  }
1186  }
1187 
1188  /// \short The number of pressure nodes is 2^DIM. The number of
1189  /// velocity nodes is the same as the number of geometric nodes.
1190  unsigned ninterpolating_node(const int &value_id)
1191  {
1192  //Find the number of interpolated values in the BASIC equations
1193  int n_basic_values = BASIC::ncont_interpolated_values();
1194  //If the id is below this number, we call the BASIC functon
1195  if(value_id < n_basic_values)
1196  {
1197  return BASIC::ninterpolating_node(value_id);
1198  }
1199  //Otherwise it's the solid and its value_id is the the current
1200  //it minus n_basic_values
1201  else
1202  {
1203  return SOLID::ninterpolating_node((value_id-n_basic_values));
1204  }
1205  }
1206 
1207  /// \short The basis interpolating the pressure is given by pshape().
1208  //// The basis interpolating the velocity is shape().
1210  Shape &psi,
1211  const int &value_id) const
1212  {
1213  //Find the number of interpolated values in the BASIC equations
1214  int n_basic_values = BASIC::ncont_interpolated_values();
1215  //If the id is below this number, we call the BASIC functon
1216  if(value_id < n_basic_values)
1217  {
1218  return BASIC::interpolating_basis(s,psi,value_id);
1219  }
1220  //Otherwise it's the solid and its value_id is the the current
1221  //it minus n_basic_values
1222  else
1223  {
1224  return SOLID::interpolating_basis(s,psi,(value_id-n_basic_values));
1225  }
1226  }
1227 
1228 
1229  /// \short Number of 'flux' terms for Z2 error estimation: Error estimation
1230  /// is based on error in BASIC element
1231  unsigned num_Z2_flux_terms() {return BASIC::num_Z2_flux_terms();}
1232 
1233  /// 'Flux' vector for Z2 error estimation: Error estimation
1234  /// is based on error in BASIC element
1236  {
1237  BASIC::get_Z2_flux(s,flux);
1238  }
1239 
1240  /// \short Perform additional hanging node procedures for variables
1241  /// that are not interpolated by all nodes. Done for both of the
1242  /// underlying element types.
1244  {
1245  BASIC::further_setup_hanging_nodes();
1246  SOLID::further_setup_hanging_nodes();
1247  }
1248 
1249 
1250  /// \short Build function: Call the one for the SOLID element since it
1251  /// calls the one basic build function automatically.
1252  void build(Mesh*& mesh_pt, Vector<Node*>& new_node_pt,
1253  bool& was_already_built,
1254  std::ofstream& new_nodes_file)
1255  {
1256  SOLID::build(mesh_pt,new_node_pt,
1257  was_already_built,new_nodes_file);
1258  }
1259 
1260 
1261  /// \short Build function: Call the one for the SOLID element since it
1262  /// calls the one basic build function automatically.
1263  void build(Mesh*& mesh_pt, Vector<Node*>& new_node_pt,
1264  bool& was_already_built)
1265  {
1266  SOLID::build(mesh_pt,new_node_pt,was_already_built);
1267  }
1268 
1269  /// Further build: Done for both of the
1270  /// underlying element types.
1272  {
1273  BASIC::further_build();
1274  SOLID::further_build();
1275  }
1276 
1277 
1278  /// \short Number of vertex nodes in the element
1279  unsigned nvertex_node() const
1280  {return BASIC::nvertex_node();}
1281 
1282  /// \short Pointer to the j-th vertex node in the element
1283  Node* vertex_node_pt(const unsigned& j) const
1284  {return BASIC::vertex_node_pt(j);}
1285 
1286  /// \short Plot the error when compared against a given exact flux.
1287  /// Also calculates the norm of the error and that of the exact flux.
1288  /// Use version in BASIC element
1290  std::ostream &outfile,
1292  double& error, double& norm)
1293  {
1294  BASIC::compute_exact_Z2_error(outfile, exact_flux_pt,
1295  error, norm);
1296  }
1297 
1298  /// \short Order of recovery shape functions for Z2 error estimation: Done
1299  /// for BASIC element since it determines the refinement
1300  unsigned nrecovery_order() {return BASIC::nrecovery_order();}
1301 
1302  ///Overload the output function: Use that of the BASIC element
1303  void output(std::ostream &outfile) {BASIC::output(outfile);}
1304 
1305  ///Output function, plotting at n_p points: Use that of the BASIC element
1306  void output(std::ostream &outfile, const unsigned &n_p)
1307  {BASIC::output(outfile,n_p);}
1308 
1309  ///Overload the output function: Use that of the BASIC element
1310  void output(FILE* file_pt) {BASIC::output(file_pt);}
1311 
1312  ///Output function: Use that of the BASIC element
1313  void output(FILE* file_pt, const unsigned &n_p)
1314  {BASIC::output(file_pt,n_p);}
1315 
1316  /// \short The number of "DOF types" that degrees of freedom in this element
1317  /// are sub-divided into.
1318  unsigned ndof_types() const
1319  {
1320  return BASIC::ndof_types() + SOLID::ndof_types();
1321  }
1322 
1323  /// \short return the number of DOF types associated with the BASIC
1324  /// elements in this combined element
1325  unsigned nbasic_dof_types() const
1326  {
1327  return BASIC::ndof_types();
1328  }
1329 
1330  /// \short return the number of DOF types associated with the SOLID
1331  /// elements in this combined element
1332  unsigned nsolid_dof_types() const
1333  {
1334  return SOLID::ndof_types();
1335  }
1336 
1337  /// \short Create a list of pairs for all unknowns in this element,
1338  /// so that the first entry in each pair contains the global equation
1339  /// number of the unknown, while the second one contains the number
1340  /// of the "DOF type" that this unknown is associated with.
1341  /// This method combines the get_dof_numbers_for_unknowns(...)
1342  /// method for the BASIC and SOLID elements. The basic elements
1343  /// retain their DOF type numbering and the SOLID elements
1344  /// DOF type numbers are incremented by nbasic_dof_types().
1346  std::list<std::pair<unsigned long,unsigned> >& dof_lookup_list) const
1347  {
1348  // get the solid list
1349  std::list<std::pair<unsigned long,unsigned> > solid_list;
1350  SOLID::get_dof_numbers_for_unknowns(solid_list);
1351 
1352  // get the basic list
1353  BASIC::get_dof_numbers_for_unknowns(dof_lookup_list);
1354 
1355  // get the number of basic dof types
1356  unsigned nbasic_dof_types = BASIC::ndof_types();
1357 
1358  // add the solid lookup list to the basic lookup list
1359  // incrementing the solid dof numbers by nbasic_dof_types
1360  typedef std::list<std::pair<unsigned long,unsigned> >::iterator IT;
1361  for (IT it=solid_list.begin();
1362  it!=solid_list.end();it++)
1363  {
1364  std::pair<unsigned long,unsigned> new_pair;
1365  new_pair.first = it->first;
1366  new_pair.second = it->second + nbasic_dof_types;
1367  dof_lookup_list.push_front(new_pair);
1368  }
1369  }
1370 };
1371 
1372 
1373 ///Explicit definition of the face geometry of these elements
1374 template<class BASIC, class SOLID>
1376 public virtual FaceGeometry<SOLID>
1377 {
1378  public:
1379  /// \short Constructor calls the constructor of the SolidQElement
1380  /// (Only the Intel compiler seems to need this!)
1381  FaceGeometry() : FaceGeometry<SOLID>() {}
1382 
1383  protected:
1384 };
1385 
1386 ///Explicit definition of the face geometry of these elements
1387 template<class BASIC, class SOLID>
1389  RefineablePseudoSolidNodeUpdateElement<BASIC,SOLID> > >:
1390 public virtual FaceGeometry<FaceGeometry<SOLID> >
1391 {
1392  public:
1393 
1394  /// \short Constuctor calls the constructor of the SolidQElement
1395  /// (Only the Intel compiler seems to need this!)
1397 
1398  protected:
1399 
1400 };
1401 
1402 
1403 }
1404 
1405 #endif
unsigned required_nvalue(const unsigned &n) const
The required number of values is the sum of the two.
FaceGeometry()
Constuctor calls the constructor of the SolidQElement (Only the Intel compiler seems to need this!) ...
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Final override for residuals function: adds contributions from both underlying element types...
void output(FILE *file_pt, const unsigned &n_p)
Output function is just the same as the basic equations.
void output(FILE *file_pt)
Overload the output function: Use that of the BASIC element.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Final override for the residuals function. Contributions are added from both underlying element types...
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
cstr elem_len * i
Definition: cfortran.h:607
void get_Z2_flux(const Vector< double > &s, Vector< double > &flux)
void fill_in_shape_derivatives_by_fd(DenseMatrix< double > &jacobian)
Fill in the derivatives of the BASIC equations w.r.t. the solid position dofs.
unsigned nsolid_dof_types() const
return the number of DOF types associated with the SOLID elements in this combined element ...
void perform_auxiliary_node_update_fct()
Execute auxiliary update function (if any) – this can be used to update any nodal values following th...
Definition: nodes.h:1509
void fill_in_contribution_to_jacobian_and_mass_matrix(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix)
Final override for mass matrix function: contributions are included from both the underlying element ...
void output(std::ostream &outfile, const unsigned &n_p)
Output function: Plot at n_p plot points using the basic element's output function.
void describe_local_dofs(std::ostream &out, const std::string &current_string) const
Function to describe the local dofs of the element. The ostream specifies the output stream to which ...
void fill_in_contribution_to_jacobian_and_mass_matrix(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix)
void interpolating_basis(const Vector< double > &s, Shape &psi, const int &value_id) const
The basis interpolating the pressure is given by pshape(). / The basis interpolating the velocity is ...
void evaluate_shape_derivs_by_chain_rule()
Evaluate shape derivatives by chain rule.
char t
Definition: cfortran.h:572
unsigned ndof_types() const
The number of "DOF types" that degrees of freedom in this element are sub-divided into...
Node * interpolating_node_pt(const unsigned &n, const int &value_id)
We must compose the underlying interpolating nodes from the BASIC and SOLID equations, the BASIC ones are first.
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
unsigned required_nvalue(const unsigned &n) const
The required number of values is the sum of the two.
HangInfo *const & hanging_pt() const
Return pointer to hanging node data (this refers to the geometric hanging node status) (const version...
Definition: nodes.h:1148
unsigned nbasic_dof_types() const
return the number of DOF types associated with the BASIC elements in this combined element ...
int solid_p_nodal_index() const
We assume that the solid stuff is stored at the end of the nodes, i.e. its index is the number of con...
void assign_additional_local_eqn_numbers()
Final override for the assign__additional_local_eqn_numbers(): Call the version for the BASIC element...
void compute_exact_Z2_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_flux_pt, double &error, double &norm)
Plot the error when compared against a given exact flux. Also calculates the norm of the error and th...
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation: Done for BASIC element since it determines...
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Final override for jacobian function: Contributions are included from both the underlying element typ...
PseudoSolidNodeUpdateElement()
Constructor, call the BASIC and SOLID elements' constructors and set the "density" parameter for soli...
void get_interpolated_values(const Vector< double > &s, Vector< double > &values)
Call get_interpolated_values(...) for both of the underlying element types.
void fill_in_shape_derivatives(DenseMatrix< double > &jacobian)
Fill in the shape derivatives of the BASIC equations w.r.t. the solid position dofs.
void build(Mesh *&mesh_pt, Vector< Node * > &new_node_pt, bool &was_already_built)
Build function: Call the one for the SOLID element since it calls the one basic build function automa...
unsigned nmaster() const
Return the number of master nodes.
Definition: nodes.h:733
unsigned nsolid_dof_types() const
return the number of DOF types associated with the SOLID elements in this combined element ...
unsigned nbasic_dof_types() const
return the number of DOF types associated with the BASIC elements in this combined element ...
unsigned nvertex_node() const
Number of vertex nodes in the element.
double Zero
Static variable to hold the default value for the pseudo-solid's inertia parameter Lambda^2...
bool is_hanging() const
Test whether the node is geometrically hanging.
Definition: nodes.h:1207
void compute_exact_Z2_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_flux_pt, double &error, double &norm)
Plot the error when compared against a given exact flux. Also calculates the norm of the error and th...
A Rank 3 Tensor class.
Definition: matrices.h:1337
void output(FILE *file_pt)
Overload the output function: Call that of the basic element.
void output(std::ostream &outfile)
Overload the output function: Use that of the BASIC element.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Add the element's contribution to its residual vector and the element Jacobian matrix (wrapper) ...
void describe_local_dofs(std::ostream &out, const std::string &current_string) const
Function to describe the local dofs of the element. The ostream specifies the output stream to which ...
void build(Mesh *&mesh_pt, Vector< Node * > &new_node_pt, bool &was_already_built, std::ofstream &new_nodes_file)
Build function: Call the one for the SOLID element since it calls the one basic build function automa...
FaceGeometry()
Constuctor calls the constructor of the SolidQElement (Only the Intel compiler seems to need this!) ...
void output(std::ostream &outfile)
Overload the output function: Call that of the basic element.
static char t char * s
Definition: cfortran.h:572
void get_dof_numbers_for_unknowns(std::list< std::pair< unsigned long, unsigned > > &dof_lookup_list) const
Create a list of pairs for all unknowns in this element, so that the first entry in each pair contain...
void rebuild_from_sons(Mesh *&mesh_pt)
Call rebuild_from_sons() for both of the underlying element types.
unsigned ndof_types() const
The number of "DOF types" that degrees of freedom in this element are sub-divided into...
void get_interpolated_values(const unsigned &t, const Vector< double > &s, Vector< double > &values)
Call get_interpolated_values(...) for both of the underlying element types.
int solid_p_nodal_index() const
We assume that the solid stuff is stored at the end of the nodes, i.e. its index is the number of con...
double local_one_d_fraction_of_interpolating_node(const unsigned &n1d, const unsigned &i, const int &value_id)
The pressure nodes are the corner nodes, so when value_id==0, the fraction is the same as the 1d node...
void output(FILE *file_pt, const unsigned &n_p)
Output function: Use that of the BASIC element.
unsigned ncont_interpolated_values() const
The number of continuously interpolated values is the sum of the SOLID and BASIC values.
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as ...
Definition: elements.h:1720
Node * get_interpolating_node_at_local_coordinate(const Vector< double > &s, const int &value_id)
The velocity nodes are the same as the geometric nodes. The pressure nodes must be calculated by usin...
double & x_gen(const unsigned &k, const unsigned &i)
Reference to the generalised position x(k,i).
Definition: nodes.h:1055
FaceGeometry()
Constructor calls the constructor of the SolidQElement (Only the Intel compiler seems to need this!) ...
void get_Z2_flux(const Vector< double > &s, Vector< double > &flux)
void output(std::ostream &outfile)
Output with default number of plot points.
void fill_in_shape_derivatives_by_fd(DenseMatrix< double > &jacobian)
Fill in the derivatives of the BASIC equations w.r.t. to the solid position dofs, taking hanging node...
unsigned num_Z2_flux_terms()
Number of 'flux' terms for Z2 error estimation: Error estimation is based on error in BASIC element...
Class that contains data for hanging nodes.
Definition: nodes.h:684
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Add the element's contribution to its residual vector (wrapper)
Refineable version of the PseudoSolidNodeUpdateELement.
unsigned ninterpolating_node(const int &value_id)
The number of pressure nodes is 2^DIM. The number of velocity nodes is the same as the number of geom...
void get_dof_numbers_for_unknowns(std::list< std::pair< unsigned long, unsigned > > &dof_lookup_list) const
Create a list of pairs for all unknowns in this element, so that the first entry in each pair contain...
void identify_geometric_data(std::set< Data * > &geometric_data_pt)
Specify Data that affects the geometry of the element by adding the position Data to the set that's p...
FaceGeometry()
Constuctor calls the constructor of the SolidQElement (Only the Intel compiler seems to need this!) ...
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation: Done for BASIC element since it determines...
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Final override for jacobian function: Calls get_jacobian() for both of the underlying element types...
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
void output(std::ostream &outfile, const unsigned &n_p)
Output function, plotting at n_p points: Use that of the BASIC element.
Node *const & master_node_pt(const unsigned &i) const
Return a pointer to the i-th master node.
Definition: nodes.h:736
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
RefineablePseudoSolidNodeUpdateElement()
Constructor, call the BASIC and SOLID elements' constructors and set the "density" parameter for soli...
unsigned num_Z2_flux_terms()
Number of 'flux' terms for Z2 error estimation: Error estimation is based on error in BASIC element...
void compute_norm(double &el_norm)
Compute norm of solution: use the version in the BASIC class if there's any ambiguity.
void further_setup_hanging_nodes()
Perform additional hanging node procedures for variables that are not interpolated by all nodes...
void evaluate_shape_derivs_by_direct_fd()
Evaluate shape derivatives by direct finite differencing.
unsigned nvertex_node() const
Number of vertex nodes in the element.
void identify_geometric_data(std::set< Data * > &geometric_data_pt)
Specify Data that affects the geometry of the element by adding the position Data to the set that's p...
unsigned ninterpolating_node_1d(const int &value_id)
The number of 1d pressure nodes is 2, otherwise we have the positional nodes.
A general mesh class.
Definition: mesh.h:74
void fill_in_contribution_to_jacobian_and_mass_matrix(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix)
Final override for mass matrix function: contributions are included from both the underlying element ...
bool Shape_derivs_by_direct_fd
Boolean flag to indicate shape derivative method.