refineable_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 member functions for refineable elements
31 
32 //oomph-lib includes
33 #include "refineable_elements.h"
34 #include "shape.h"
35 
36 namespace oomph
37 {
38  //=====================================================================
39  /// Destructor (needed here because of IBM xlC compiler under AIX)
40  /// Delete the storage allocated for the local equations.
41  //====================================================================
43 
44 //========================================================================
45 /// Max. allowed discrepancy in element integrity check
46 //========================================================================
48 
49 //=========================================================================
50 /// Helper function that is used to check that the value_id is in the range
51 /// allowed by the element. The number of continuously interpolated values
52 /// and the value_id must be passed as arguments.
53 //========================================================================
55  &n_continuously_interpolated_values,
56  const int &value_id)
57  {
58  //If the value_id is more than the number of continuously interpolated
59  //values or less than -1 (for the position), we have a problem
60  if((value_id < -1) || (value_id >= n_continuously_interpolated_values))
61  {
62  std::ostringstream error_stream;
63  error_stream << "Value_id " << value_id << " is out of range."
64  << std::endl
65  << "It can only take the values -1 (position) "
66  << "or an integer in the range 0 to "
67  << n_continuously_interpolated_values
68  << std::endl;
69 
70  throw OomphLibError(error_stream.str(),
71  OOMPH_CURRENT_FUNCTION,
72  OOMPH_EXCEPTION_LOCATION);
73  }
74  }
75 
76 //=========================================================================
77 /// Internal function that is used to assemble the jacobian of the mapping
78 /// from local coordinates (s) to the eulerian coordinates (x), given the
79 /// derivatives of the shape functions.
80 //=========================================================================
83  DenseMatrix<double> &jacobian) const
84  {
85  //Locally cache the elemental dimension
86  const unsigned el_dim = dim();
87  //The number of shape functions must be equal to the number
88  //of nodes (by definition)
89  const unsigned n_shape = nnode();
90  //The number of shape function types must be equal to the number
91  //of nodal position types (by definition)
92  const unsigned n_shape_type = nnodal_position_type();
93 
94 #ifdef PARANOID
95  //Check for dimensional compatibility
96  if(dim() != nodal_dimension())
97  {
98  std::ostringstream error_message;
99  error_message << "Dimension mismatch" << std::endl;
100  error_message << "The elemental dimension: " << dim()
101  << " must equal the nodal dimension: "
102  << nodal_dimension()
103  << " for the jacobian of the mapping to be well-defined"
104  << std::endl;
105  throw OomphLibError(
106  error_message.str(),
107  OOMPH_CURRENT_FUNCTION,
108  OOMPH_EXCEPTION_LOCATION);
109  }
110 #endif
111 
112  //Loop over the rows of the jacobian
113  for(unsigned i=0;i<el_dim;i++)
114  {
115  //Loop over the columns of the jacobian
116  for(unsigned j=0;j<el_dim;j++)
117  {
118  //Zero the entry
119  jacobian(i,j) = 0.0;
120  //Loop over the shape functions
121  for(unsigned l=0;l<n_shape;l++)
122  {
123  for(unsigned k=0;k<n_shape_type;k++)
124  {
125  //Jacobian is dx_j/ds_i, which is represented by the sum
126  //over the dpsi/ds_i of the nodal points X j
127  //Call the Hanging version of positions
128  jacobian(i,j) += nodal_position_gen(l,k,j)*dpsids(l,k,i);
129  }
130  }
131  }
132  }
133  }
134 
135 //=========================================================================
136 /// Internal function that is used to assemble the jacobian of second
137 /// derivatives of the the mapping from local coordinates (s) to the
138 /// eulerian coordinates (x), given the second derivatives of the
139 /// shape functions.
140 //=========================================================================
143  DenseMatrix<double> &jacobian2) const
144  {
145  //Find the the dimension of the element
146  const unsigned el_dim = dim();
147  //Find the number of shape functions and shape functions types
148  //Must be equal to the number of nodes and their position types
149  //by the definition of the shape function.
150  const unsigned n_shape = nnode();
151  const unsigned n_shape_type = nnodal_position_type();
152  //Find the number of second derivatives
153  const unsigned n_row = N2deriv[el_dim];
154 
155  //Assemble the "jacobian" (d^2 x_j/ds_i^2) for second derivatives of
156  //shape functions
157  //Loop over the rows (number of second derivatives)
158  for(unsigned i=0;i<n_row;i++)
159  {
160  //Loop over the columns (element dimension
161  for(unsigned j=0;j<el_dim;j++)
162  {
163  //Zero the entry
164  jacobian2(i,j) = 0.0;
165  //Loop over the shape functions
166  for(unsigned l=0;l<n_shape;l++)
167  {
168  //Loop over the shape function types
169  for(unsigned k=0;k<n_shape_type;k++)
170  {
171  //Add the terms to the jacobian entry
172  //Call the Hanging version of positions
173  jacobian2(i,j) += nodal_position_gen(l,k,j)*d2psids(l,k,i);
174  }
175  }
176  }
177  }
178  }
179 
180  //=====================================================================
181  /// Assemble the covariant Eulerian base vectors and return them in the
182  /// matrix interpolated_G. The derivatives of the shape functions with
183  /// respect to the local coordinate should already have been calculated
184  /// before calling this function
185  //=====================================================================
187  const DShape &dpsids, DenseMatrix<double> &interpolated_G) const
188  {
189  //Find the number of nodes and position types
190  const unsigned n_node = nnode();
191  const unsigned n_position_type = nnodal_position_type();
192  //Find the dimension of the node and element
193  const unsigned n_dim_node = nodal_dimension();
194  const unsigned n_dim_element = dim();
195 
196  //Loop over the dimensions and compute the entries of the
197  //base vector matrix
198  for(unsigned i=0;i<n_dim_element;i++)
199  {
200  for(unsigned j=0;j<n_dim_node;j++)
201  {
202  //Initialise the j-th component of the i-th base vector to zero
203  interpolated_G(i,j) = 0.0;
204  for(unsigned l=0;l<n_node;l++)
205  {
206  for(unsigned k=0;k<n_position_type;k++)
207  {
208  interpolated_G(i,j) += nodal_position_gen(l,k,j)*dpsids(l,k,i);
209  }
210  }
211  }
212  }
213  }
214 
215 
216 //==========================================================================
217 /// Calculate the mapping from local to eulerian coordinates
218 /// assuming that the coordinates are aligned in the direction of the local
219 /// coordinates, i.e. there are no cross terms and the jacobian is diagonal.
220 /// The local derivatives are passed as dpsids and the jacobian and
221 /// inverse jacobian are returned.
222 //==========================================================================
223  double RefineableElement::
225  DenseMatrix<double> &jacobian,
226  DenseMatrix<double> &inverse_jacobian)
227  const
228  {
229  //Find the dimension of the element
230  const unsigned el_dim = dim();
231  //Find the number of shape functions and shape functions types
232  //Equal to the number of nodes and their position types by definition
233  const unsigned n_shape = nnode();
234  const unsigned n_shape_type = nnodal_position_type();
235 
236 #ifdef PARANOID
237  //Check for dimension compatibility
238  if(dim() != nodal_dimension())
239  {
240  std::ostringstream error_message;
241  error_message << "Dimension mismatch" << std::endl;
242  error_message << "The elemental dimension: " << dim()
243  << " must equal the nodal dimension: "
244  << nodal_dimension()
245  << " for the jacobian of the mapping to be well-defined"
246  << std::endl;
247  throw OomphLibError(error_message.str(),
248  OOMPH_CURRENT_FUNCTION,
249  OOMPH_EXCEPTION_LOCATION);
250  }
251 #endif
252 
253  //In this case we assume that there are no cross terms, that is
254  //global coordinate 0 is always in the direction of local coordinate 0
255 
256  //Loop over the coordinates
257  for(unsigned i=0;i<el_dim;i++)
258  {
259  //Zero the jacobian and inverse jacobian entries
260  for(unsigned j=0;j<el_dim;j++)
261  {jacobian(i,j) = 0.0; inverse_jacobian(i,j) = 0.0;}
262 
263  //Loop over the shape functions
264  for(unsigned l=0;l<n_shape;l++)
265  {
266  //Loop over the types of dof
267  for(unsigned k=0;k<n_shape_type;k++)
268  {
269  //Derivatives are always dx_{i}/ds_{i}
270  jacobian(i,i) += nodal_position_gen(l,k,i)*dpsids(l,k,i);
271  }
272  }
273  }
274 
275  //Now calculate the determinant of the matrix
276  double det = 1.0;
277  for(unsigned i=0;i<el_dim;i++) {det *= jacobian(i,i);}
278 
279 //Report if Matrix is singular, or negative
280 #ifdef PARANOID
281  check_jacobian(det);
282 #endif
283 
284  //Calculate the inverse mapping (trivial in this case)
285  for(unsigned i=0;i<el_dim;i++)
286  {inverse_jacobian(i,i) = 1.0/jacobian(i,i);}
287 
288  //Return the value of the Jacobian
289  return(det);
290  }
291 
292 //========================================================================
293 /// Deactivate the element by marking setting all local pointers to
294 /// obsolete nodes to zero
295 //=======================================================================
297  {
298  //Find the number of nodes
299  const unsigned n_node = nnode();
300  //Loop over the nodes
301  for(unsigned n=0;n<n_node;n++)
302  {
303  //If the node pointer has not already been nulled, but is obsolete
304  //then null it
305  if((this->node_pt(n)!=0) && (this->node_pt(n)->is_obsolete()))
306  {this->node_pt(n) = 0;}
307  }
308  }
309 
310 //=======================================================================
311 /// Assign the local hang eqn.
312 //=======================================================================
314  const bool &store_local_dof_pt)
315  {
316  //Find the number of nodes
317  const unsigned n_node = nnode();
318 
319  //Check there are nodes!
320  if(n_node > 0)
321  {
322  //Find the number of continuously interpolated values
323  const unsigned n_cont_values = ncont_interpolated_values();
324 
325  //Delete the storage associated with the local hanging equation schemes
326  delete[] Local_hang_eqn;
327  //Allocate new storage
328  Local_hang_eqn = new std::map<Node*,int>[n_cont_values];
329 
330  //Boolean that is set to true if there are hanging equation numbers
331  bool hanging_eqn_numbers = false;
332 
333  //Maps that store whether the node's local equation number for a
334  //particular value has already been assigned
335  Vector<std::map<Node*,bool> > local_eqn_number_done(n_cont_values);
336 
337  //Get number of dofs assigned thus far
338  unsigned local_eqn_number = ndof();
339 
340  //A local queue to store the global equation numbers
341  std::deque<unsigned long> global_eqn_number_queue;
342 
343  //Now loop over all the nodes again to find the master nodes
344  //external to the element
345  for(unsigned n=0;n<n_node;n++)
346  {
347  //Loop over the number of continuously-interpolated values at the node
348  for(unsigned j=0;j<n_cont_values;j++)
349  {
350  //If the node is hanging node in value j
351  if(node_pt(n)->is_hanging(j))
352  {
353  //Get the pointer to the appropriate hanging info object
354  HangInfo* hang_info_pt = node_pt(n)->hanging_pt(j);
355  //Find the number of master nodes
356  unsigned n_master = hang_info_pt->nmaster();
357  //Loop over the master nodes
358  for(unsigned m=0;m<n_master;m++)
359  {
360  //Get the m-th master node
361  Node* Master_node_pt = hang_info_pt->master_node_pt(m);
362 
363  //If the master node's value has not been considered already,
364  //give it a local equation number
365  if(local_eqn_number_done[j][Master_node_pt] == false)
366  {
367 #ifdef PARANOID
368  //Check that the value is stored at the master node
369  //If not the master node must have be set up incorrectly
370  if(Master_node_pt->nvalue() < j)
371  {
372  std::ostringstream error_stream;
373  error_stream
374  << "Master node for " << j
375  << "-th value only has " << Master_node_pt->nvalue()
376  << " stored values!" << std::endl;
377 
378  throw OomphLibError(
379  error_stream.str(),
380  OOMPH_CURRENT_FUNCTION,
381  OOMPH_EXCEPTION_LOCATION);
382  }
383 #endif
384 
385  //We now need to test whether the master node is actually
386  //a local node, in which case its local equation number
387  //will already have been assigned and will be stored
388  //in nodal_local_eqn
389 
390  //Storage for the index of the local node
391  //Initialised to n_node (beyond the possible indices of
392  //"real" nodes)
393  unsigned local_node_index=n_node;
394  //Loop over the local nodes (again)
395  for(unsigned n1=0;n1<n_node;n1++)
396  {
397  //If the master node is a local node
398  //get its index and break out of the loop
399  if(Master_node_pt==node_pt(n1))
400  {
401  local_node_index=n1;
402  break;
403  }
404  }
405 
406  //Now we test whether the node was found
407  if(local_node_index < n_node)
408  {
409  //Copy the local equation number to the
410  //pointer-based look-up scheme
411  Local_hang_eqn[j][Master_node_pt] =
412  nodal_local_eqn(local_node_index,j);
413  }
414  //Otherwise it's a new master node
415  else
416  {
417  //Get the equation number
418  long eqn_number = Master_node_pt->eqn_number(j);
419  //If equation_number positive add to array
420  if(eqn_number >= 0)
421  {
422  //Add global equation number to the queue
423  global_eqn_number_queue.push_back(eqn_number);
424  //Add pointer to the dof to the queue if required
425  if(store_local_dof_pt)
426  {
428  Master_node_pt->value_pt(j));
429  }
430  //Add to pointer based scheme
431  Local_hang_eqn[j][Master_node_pt] = local_eqn_number;
432  //Increase number of local variables
433  local_eqn_number++;
434  }
435  //Otherwise the value is pinned
436  else
437  {
438  Local_hang_eqn[j][Master_node_pt] = Data::Is_pinned;
439  }
440  //There are now hanging equation numbers
441  }
442  //The value at this node has now been done
443  local_eqn_number_done[j][Master_node_pt] = true;
444  //There are hanging equation numbers
445  hanging_eqn_numbers = true;
446  }
447  }
448  }
449  }
450  } //End of second loop over nodes
451 
452  //Now add our global equations numbers to the internal element storage
453  add_global_eqn_numbers(global_eqn_number_queue,
455  //Clear the memory used in the deque
456  if(store_local_dof_pt)
457  {std::deque<double*>().swap(GeneralisedElement::Dof_pt_deque);}
458 
459 
460  //If there are no hanging_eqn_numbers delete the (empty) stored maps
461  if(!hanging_eqn_numbers)
462  {
463  delete[] Local_hang_eqn;
464  Local_hang_eqn=0;
465  }
466 
467 
468  // Setup map that associates a unique number with any of the nodes
469  // that actively control the shape of the element (i.e. they are
470  // either non-hanging nodes of this element or master nodes
471  // of hanging nodes.
472  unsigned count=0;
473  std::set<Node*> all_nodes;
475  for (unsigned j=0;j<n_node;j++)
476  {
477  // Get node
478  Node* nod_pt=node_pt(j);
479 
480  //If the node is hanging, consider master nodes
481  if(nod_pt->is_hanging())
482  {
483  HangInfo* hang_info_pt = node_pt(j)->hanging_pt();
484  unsigned n_master = hang_info_pt->nmaster();
485 
486  //Loop over the master nodes
487  for(unsigned m=0;m<n_master;m++)
488  {
489  Node* master_node_pt=hang_info_pt->master_node_pt(m);
490  // Do we have this one already?
491  unsigned old_size=all_nodes.size();
492  all_nodes.insert(master_node_pt);
493  if (all_nodes.size()>old_size)
494  {
495  Shape_controlling_node_lookup[master_node_pt]=count;
496  count++;
497  }
498  }
499  }
500  // Not hanging: Consider the node itself
501  else
502  {
503  // Do we have this one already?
504  unsigned old_size=all_nodes.size();
505  all_nodes.insert(nod_pt);
506  if (all_nodes.size()>old_size)
507  {
508  Shape_controlling_node_lookup[nod_pt]=count;
509  count++;
510  }
511  }
512  }
513 
514  } //End of if nodes
515  }
516 
517  //======================================================================
518  /// The purpose of this function is to identify all possible
519  /// Data that can affect the fields interpolated by the FiniteElement.
520  /// This must be overloaded to include data from any hanging nodes
521  /// correctly
522  //=======================================================================
524  std::set<std::pair<Data*,unsigned> > &paired_field_data)
525 {
526  //Loop over all internal data
527  const unsigned n_internal = this->ninternal_data();
528  for(unsigned n=0;n<n_internal;n++)
529  {
530  //Cache the data pointer
531  Data* const dat_pt = this->internal_data_pt(n);
532  //Find the number of data values stored in the data object
533  const unsigned n_value = dat_pt->nvalue();
534  //Add the index of each data value and the pointer to the set
535  //of pairs
536  for(unsigned i=0;i<n_value;i++)
537  {
538  paired_field_data.insert(std::make_pair(dat_pt,i));
539  }
540  }
541 
542  //Loop over all the nodes
543  const unsigned n_node = this->nnode();
544  for(unsigned n=0;n<n_node;n++)
545  {
546  //Cache the node pointer
547  Node* const nod_pt = this->node_pt(n);
548  //Find the number of values stored at the node
549  const unsigned n_value = nod_pt->nvalue();
550 
551  //Loop over the number of values
552  for(unsigned i=0;i<n_value;i++)
553  {
554  //If the node is hanging in the variable
555  if(nod_pt->is_hanging(i))
556  {
557  //Cache the pointer to the HangInfo object for this variable
558  HangInfo* const hang_pt = nod_pt->hanging_pt(i);
559  //Get number of master nodes
560  const unsigned nmaster = hang_pt->nmaster();
561 
562  // Loop over masters
563  for(unsigned m=0;m<nmaster;m++)
564  {
565  //Cache the pointer to the master node
566  Node* const master_nod_pt=hang_pt->master_node_pt(m);
567 
568  //Under the assumption that the same data value is interpolated
569  //by the hanging nodes, which it really must be
570  //Add it to the paired field data
571  paired_field_data.insert(std::make_pair(master_nod_pt,i));
572  }
573  }
574  //Otherwise the node is not hanging in the variable, treat as normal
575  else
576  {
577  //Add the index of each data value and the pointer to the set
578  //of pairs
579  paired_field_data.insert(std::make_pair(nod_pt,i));
580  }
581  }
582  }
583 }
584 
585 
586 
587 //==========================================================================
588 /// Compute derivatives of elemental residual vector with respect
589 /// to nodal coordinates. Default implementation by FD can be overwritten
590 /// for specific elements.
591 /// dresidual_dnodal_coordinates(l,i,j) = d res(l) / dX_{ij}
592 /// This version is overloaded from the version in FiniteElement
593 /// and takes hanging nodes into account -- j in the above loop
594 /// loops over all the nodes that actively control the
595 /// shape of the element (i.e. they are non-hanging or master nodes of
596 /// hanging nodes in this element).
597 //==========================================================================
599  RankThreeTensor<double>& dresidual_dnodal_coordinates)
600  {
601  // Number of nodes
602  unsigned n_nod=nnode();
603 
604  // If the element has no nodes (why??!!) return straightaway
605  if (n_nod==0) return;
606 
607  // Get dimension from first node
608  unsigned dim_nod=node_pt(0)->ndim();
609 
610  // Number of dofs
611  unsigned n_dof=ndof();
612 
613  // Get reference residual
614  Vector<double> res(n_dof);
615  Vector<double> res_pls(n_dof);
616  get_residuals(res);
617 
618  // FD step
620 
621  // Do FD loop over all active nodes
622  for (std::map<Node*,unsigned>::iterator it=
625  it++)
626  {
627  // Get node
628  Node* nod_pt=it->first;
629 
630  // Get its number
631  unsigned node_number=it->second;
632 
633  // Loop over coordinate directions
634  for (unsigned i=0;i<dim_nod;i++)
635  {
636  // Make backup
637  double backup=nod_pt->x(i);
638 
639  // Do FD step. No node update required as we're
640  // attacking the coordinate directly...
641  nod_pt->x(i)+=eps_fd;
642 
643  // Perform auxiliary node update function
645 
646  // Get advanced residual
647  get_residuals(res_pls);
648 
649  // Fill in FD entries [Loop order is "wrong" here as l is the
650  // slow index but this is in a function that's costly anyway
651  // and gives us the fastest loop outside where these tensor
652  // is actually used.]
653  for (unsigned l=0;l<n_dof;l++)
654  {
655  dresidual_dnodal_coordinates(l,i,node_number)=
656  (res_pls[l]-res[l])/eps_fd;
657  }
658 
659  // Reset coordinate. No node update required as we're
660  // attacking the coordinate directly...
661  nod_pt->x(i)=backup;
662 
663  // Perform auxiliary node update function
665 
666  }
667  }
668 
669  }
670 
671 
672 //============================================================================
673 /// This function calculates the entries of Jacobian matrix, used in
674 /// the Newton method, associated with the nodal degrees of freedom.
675 /// This is done by finite differences to handle the general case
676 /// Overload the standard case to include hanging node case
677 //==========================================================================
680  DenseMatrix<double> &jacobian)
681  {
682  //Find the number of nodes
683  const unsigned n_node = nnode();
684 
685  //If there are no nodes, return straight away
686  if(n_node == 0) {return;}
687 
688  //Call the update function to ensure that the element is in
689  //a consistent state before finite differencing starts
691 
692 // bool use_first_order_fd=false;
693 
694  //Find the number of dofs in the element
695  const unsigned n_dof = ndof();
696 
697  //Create newres vector
698  Vector<double> newres(n_dof); // , newres_minus(n_dof);
699 
700  //Used default value defined in GeneralisedElement
701  const double fd_step = Default_fd_jacobian_step;
702 
703  //Integer storage for local unknowns
704  int local_unknown=0;
705 
706  //Loop over the nodes
707  for(unsigned n=0;n<n_node;n++)
708  {
709  //Get the number of values stored at the node
710  const unsigned n_value = node_pt(n)->nvalue();
711 
712  //Get a pointer to the local node
713  Node* const local_node_pt = node_pt(n);
714 
715  //Loop over the number of values stored at the node
716  for(unsigned i=0;i<n_value;i++)
717  {
718  //If the nodal value is not hanging
719  if(local_node_pt->is_hanging(i)==false)
720  {
721  //Get the equation number
722  local_unknown = nodal_local_eqn(n,i);
723  //If the variable is free
724  if(local_unknown >= 0)
725  {
726  //Get a pointer to the nodal value
727  double* const value_pt = local_node_pt->value_pt(i);
728 
729  //Save the old value of nodal value
730  const double old_var = *value_pt;
731 
732  //Increment the nodal value
733  *value_pt += fd_step;
734 
735  //Now update any slaved variables
737 
738  //Calculate the new residuals
739  get_residuals(newres);
740 
741 // if (use_first_order_fd)
742  {
743  //Do forward finite differences
744  for(unsigned m=0;m<n_dof;m++)
745  {
746  //Stick the entry into the Jacobian matrix
747  jacobian(m,local_unknown) = (newres[m] - residuals[m])/fd_step;
748  }
749  }
750 // else
751 // {
752 // //Take backwards step
753 // *value_pt = old_var - fd_step;
754 
755 // //Calculate the new residuals at backward position
756 // get_residuals(newres_minus);
757 
758 // //Do central finite differences
759 // for(unsigned m=0;m<n_dof;m++)
760 // {
761 // //Stick the entry into the Jacobian matrix
762 // jacobian(m,local_unknown) =
763 // (newres[m] - newres_minus[m])/(2.0*fd_step);
764 // }
765 // }
766 
767  //Reset the nodal value
768  *value_pt = old_var;
769 
770  //Reset any slaved variables
772  }
773  }
774  //Otherwise the value is hanging node
775  else
776  {
777  //Get the local hanging infor
778  HangInfo *hang_info_pt = local_node_pt->hanging_pt(i);
779  //Loop over the master nodes
780  const unsigned n_master = hang_info_pt->nmaster();
781  for(unsigned m=0;m<n_master;m++)
782  {
783  //Get the pointer to the master node
784  Node* const master_node_pt = hang_info_pt->master_node_pt(m);
785 
786  //Get the number of the unknown
787  local_unknown = local_hang_eqn(master_node_pt,i);
788  //If the variable is free
789  if(local_unknown >= 0)
790  {
791  //Get a pointer to the nodal value stored at the master node
792  double* const value_pt = master_node_pt->value_pt(i);
793 
794  //Save the old value of the nodal value stored at the master node
795  const double old_var = *value_pt;
796 
797  //Increment the nodal value stored at the master node
798  *value_pt += fd_step;
799 
800  //Now update any slaved variables
802 
803  //Calculate the new residuals
804  get_residuals(newres);
805 
806 // if (use_first_order_fd)
807  {
808  //Do forward finite differences
809  for(unsigned m=0;m<n_dof;m++)
810  {
811  //Stick the entry into the Jacobian matrix
812  jacobian(m,local_unknown) = (newres[m] - residuals[m])/fd_step;
813  }
814  }
815 // else
816 // {
817 // //Take backwards step
818 // *value_pt = old_var - fd_step;
819 
820 // //Calculate the new residuals at backward position
821 // get_residuals(newres_minus);
822 
823 // //Do central finite differences
824 // for(unsigned m=0;m<n_dof;m++)
825 // {
826 // //Stick the entry into the Jacobian matrix
827 // jacobian(m,local_unknown) =
828 // (newres[m] - newres_minus[m])/(2.0*fd_step);
829 // }
830 // }
831 
832  //Reset the value at the master node
833  *value_pt = old_var;
834 
835  //Reset any slaved variables
837  }
838  } //End of loop over master nodes
839  } //End of hanging node case
840  } //End of loop over values
841  } //End of loop over nodes
842 
843  //End of finite difference loop
844  //Final reset of any slaved data
846  }
847 
848 
849 //=========================================================================
850 /// Internal function that is used to assemble the jacobian of the mapping
851 /// from local coordinates (s) to the lagrangian coordinates (xi), given the
852 /// derivatives of the shape functions.
853 //=========================================================================
856  DenseMatrix<double> &jacobian) const
857  {
858  //Find the the dimension of the element
859  const unsigned el_dim = dim();
860  //Find the number of shape functions and shape functions types
861  const unsigned n_shape = nnode();
862  const unsigned n_shape_type = nnodal_lagrangian_type();
863 
864  //Loop over the rows of the jacobian
865  for(unsigned i=0;i<el_dim;i++)
866  {
867  //Loop over the columns of the jacobian
868  for(unsigned j=0;j<el_dim;j++)
869  {
870  //Zero the entry
871  jacobian(i,j) = 0.0;
872  //Loop over the shape functions
873  for(unsigned l=0;l<n_shape;l++)
874  {
875  for(unsigned k=0;k<n_shape_type;k++)
876  {
877  //Jacobian is dx_j/ds_i, which is represented by the sum
878  //over the dpsi/ds_i of the nodal points X j
879  //Use the hanging version here
880  jacobian(i,j) += lagrangian_position_gen(l,k,j)*dpsids(l,k,i);
881  }
882  }
883  }
884  }
885  }
886 
887 //=========================================================================
888 /// Internal function that is used to assemble the jacobian of second
889 /// derivatives of the the mapping from local coordinates (s) to the
890 /// lagrangian coordinates (xi), given the second derivatives of the
891 /// shape functions.
892 //=========================================================================
895  DenseMatrix<double> &jacobian2) const
896  {
897  //Find the the dimension of the element
898  const unsigned el_dim = dim();
899  //Find the number of shape functions and shape functions types
900  const unsigned n_shape = nnode();
901  const unsigned n_shape_type = nnodal_lagrangian_type();
902  //Find the number of second derivatives
903  const unsigned n_row = N2deriv[el_dim];
904 
905  //Assemble the "jacobian" (d^2 x_j/ds_i^2) for second derivatives of
906  //shape functions
907  //Loop over the rows (number of second derivatives)
908  for(unsigned i=0;i<n_row;i++)
909  {
910  //Loop over the columns (element dimension
911  for(unsigned j=0;j<el_dim;j++)
912  {
913  //Zero the entry
914  jacobian2(i,j) = 0.0;
915  //Loop over the shape functions
916  for(unsigned l=0;l<n_shape;l++)
917  {
918  //Loop over the shape function types
919  for(unsigned k=0;k<n_shape_type;k++)
920  {
921  //Add the terms to the jacobian entry
922  //Use the hanging version here
923  jacobian2(i,j) += lagrangian_position_gen(l,k,j)*d2psids(l,k,i);
924  }
925  }
926  }
927  }
928  }
929 
930 //==========================================================================
931 /// Calculate the mapping from local to lagrangian coordinates
932 /// assuming that the coordinates are aligned in the direction of the local
933 /// coordinates, i.e. there are no cross terms and the jacobian is diagonal.
934 /// The local derivatives are passed as dpsids and the jacobian and
935 /// inverse jacobian are returned.
936 //==========================================================================
939  DenseMatrix<double> &jacobian,
940  DenseMatrix<double> &inverse_jacobian)
941  const
942  {
943  //Find the dimension of the element
944  const unsigned el_dim = dim();
945  //Find the number of shape functions and shape functions types
946  const unsigned n_shape = nnode();
947  const unsigned n_shape_type = nnodal_lagrangian_type();
948 
949  //In this case we assume that there are no cross terms, that is
950  //global coordinate 0 is always in the direction of local coordinate 0
951 
952  //Loop over the coordinates
953  for(unsigned i=0;i<el_dim;i++)
954  {
955  //Zero the jacobian and inverse jacobian entries
956  for(unsigned j=0;j<el_dim;j++)
957  {jacobian(i,j) = 0.0; inverse_jacobian(i,j) = 0.0;}
958 
959  //Loop over the shape functions
960  for(unsigned l=0;l<n_shape;l++)
961  {
962  //Loop over the types of dof
963  for(unsigned k=0;k<n_shape_type;k++)
964  {
965  //Derivatives are always dx_{i}/ds_{i}
966  jacobian(i,i) += lagrangian_position_gen(l,k,i)*dpsids(l,k,i);
967  }
968  }
969  }
970 
971  //Now calculate the determinant of the matrix
972  double det = 1.0;
973  for(unsigned i=0;i<el_dim;i++) {det *= jacobian(i,i);}
974 
975 //Report if Matrix is singular, or negative
976 #ifdef PARANOID
977  check_jacobian(det);
978 #endif
979 
980  //Calculate the inverse mapping (trivial in this case)
981  for(unsigned i=0;i<el_dim;i++)
982  {inverse_jacobian(i,i) = 1.0/jacobian(i,i);}
983 
984  //Return the value of the Jacobian
985  return(det);
986  }
987 
988 
989 
990 //========================================================================
991 /// The number of geometric data affecting a
992 /// RefineableSolidFiniteElement is the positional Data of all
993 /// non-hanging nodes plus the geometric Data of all distinct
994 /// master nodes. Recomputed on the fly.
995 //========================================================================
997  {
998 
999  //Find the number of nodes
1000  const unsigned n_node = nnode();
1001 
1002  // Temporary storage for unique position data
1003  std::set<Data*> all_position_data_pt;
1004 
1005  //Now loop over all the nodes again to find the master nodes
1006  //of any hanging nodes that have not yet been assigned
1007  for(unsigned n=0;n<n_node;n++)
1008  {
1009 
1010  //If the node is a hanging node
1011  if(node_pt(n)->is_hanging())
1012  {
1013  //Find the local hang info object
1014  HangInfo* hang_info_pt = node_pt(n)->hanging_pt();
1015 
1016  //Find the number of master nodes
1017  unsigned n_master = hang_info_pt->nmaster();
1018 
1019  //Loop over the master nodes
1020  for(unsigned m=0;m<n_master;m++)
1021  {
1022  //Get the m-th master node
1023  Node* Master_node_pt = hang_info_pt->master_node_pt(m);
1024 
1025  // Add to set
1026  all_position_data_pt.insert(
1027  dynamic_cast<SolidNode*>(Master_node_pt)->variable_position_pt());
1028  }
1029  }
1030  // Not hanging
1031  else
1032  {
1033  // Add node itself to set
1034  all_position_data_pt.insert(
1035  dynamic_cast<SolidNode*>(node_pt(n))->variable_position_pt());
1036  }
1037 
1038  } //End of loop over nodes
1039 
1040  // How many are there?
1041  return all_position_data_pt.size();
1042 
1043  }
1044 
1045 
1046 
1047 //========================================================================
1048 /// \short Return pointer to the j-th Data item that the object's
1049 /// shape depends on: Positional data of non-hanging nodes and
1050 /// positional data of master nodes. Recomputed on the fly.
1051 //========================================================================
1053  {
1054 
1055  //Find the number of nodes
1056  const unsigned n_node = nnode();
1057 
1058  // Temporary storage for unique position data. Set and vector are
1059  // required to ensure uniqueness in enumeration on different processors.
1060  // Set checks uniqueness; Vector stores entries in predictable order.
1061  std::set<Data*> all_position_data_pt;
1062  Vector<Data*> all_position_data_vector_pt;
1063 
1064  // Number of entries in set before possibly adding new entry
1065  unsigned n_old=0;
1066 
1067  //Now loop over all the nodes again to find the master nodes
1068  //of any hanging nodes that have not yet been assigned
1069  for(unsigned n=0;n<n_node;n++)
1070  {
1071  //If the node is a hanging node
1072  if(node_pt(n)->is_hanging())
1073  {
1074  //Find the local hang info object
1075  HangInfo* hang_info_pt = node_pt(n)->hanging_pt();
1076 
1077  //Find the number of master nodes
1078  unsigned n_master = hang_info_pt->nmaster();
1079 
1080  //Loop over the master nodes
1081  for(unsigned m=0;m<n_master;m++)
1082  {
1083  //Get the m-th master node
1084  Node* Master_node_pt = hang_info_pt->master_node_pt(m);
1085 
1086  // Positional data
1087  Data* pos_data_pt=
1088  dynamic_cast<SolidNode*>(Master_node_pt)->variable_position_pt();
1089 
1090  // Add to set
1091  n_old=all_position_data_pt.size();
1092  all_position_data_pt.insert(pos_data_pt);
1093 
1094  // New entry?
1095  if (all_position_data_pt.size()>n_old)
1096  {
1097  all_position_data_vector_pt.push_back(pos_data_pt);
1098  }
1099  }
1100  }
1101  // Not hanging
1102  else
1103  {
1104  // Add node itself to set
1105 
1106  // Positional data
1107  Data* pos_data_pt=
1108  dynamic_cast<SolidNode*>(node_pt(n))->variable_position_pt();
1109 
1110  // Add to set
1111  n_old=all_position_data_pt.size();
1112  all_position_data_pt.insert(pos_data_pt);
1113 
1114  // New entry?
1115  if (all_position_data_pt.size()>n_old)
1116  {
1117  all_position_data_vector_pt.push_back(pos_data_pt);
1118  }
1119  }
1120 
1121  } //End of loop over nodes
1122 
1123 
1124  // Return j-th entry
1125  return all_position_data_vector_pt[j];
1126 
1127  }
1128 
1129 
1130 //========================================================================
1131 /// \short Specify Data that affects the geometry of the element
1132 /// by adding the position Data to the set that's passed in.
1133 /// (This functionality is required in FSI problems; set is used to
1134 /// avoid double counting). Refineable version includes hanging nodes
1135 //========================================================================
1137  geometric_data_pt)
1138  {
1139  //Loop over the node update data and add to the set
1140  const unsigned n_node=this->nnode();
1141  for(unsigned j=0;j<n_node;j++)
1142  {
1143 
1144  //If the node is a hanging node
1145  if(node_pt(j)->is_hanging())
1146  {
1147  //Find the local hang info object
1148  HangInfo* hang_info_pt = node_pt(j)->hanging_pt();
1149 
1150  //Find the number of master nodes
1151  unsigned n_master = hang_info_pt->nmaster();
1152 
1153  //Loop over the master nodes
1154  for(unsigned m=0;m<n_master;m++)
1155  {
1156  //Get the m-th master node
1157  Node* Master_node_pt = hang_info_pt->master_node_pt(m);
1158 
1159  // Add to set
1160  geometric_data_pt.insert(
1161  dynamic_cast<SolidNode*>(Master_node_pt)->variable_position_pt());
1162  }
1163  }
1164  // Not hanging
1165  else
1166  {
1167  // Add node itself to set
1168  geometric_data_pt.insert(
1169  dynamic_cast<SolidNode*>(node_pt(j))->variable_position_pt());
1170  }
1171  }
1172  }
1173 
1174 //========================================================================
1175 /// The standard equation numbering scheme for solid positions,
1176 /// so that hanging Node information is included.
1177 //========================================================================
1179  const bool &store_local_dof_pt)
1180  {
1181  //Find the number of nodes
1182  const unsigned n_node = nnode();
1183 
1184  //Check there are nodes!
1185  if(n_node > 0)
1186  {
1187  //Wipe the local matrix maps, they will be assigned below
1188  Local_position_hang_eqn.clear();
1189 
1190  //Find the local numbers
1191  const unsigned n_position_type = nnodal_position_type();
1192  const unsigned nodal_dim = nodal_dimension();
1193 
1194  //Matrix structure to store all positional equations at a node
1195  DenseMatrix<int> Position_local_eqn_at_node(n_position_type,
1196  nodal_dim);
1197 
1198  //Map that store whether the node's equation numbers have already been
1199  //added to the local arrays
1200  std::map<Node*, bool> local_eqn_number_done;
1201 
1202  //Get number of dofs so far
1203  unsigned local_eqn_number = ndof();
1204 
1205  //A local queue to store the global equation numbers
1206  std::deque<unsigned long> global_eqn_number_queue;
1207 
1208  //Now loop over all the nodes again to find the master nodes
1209  //of any hanging nodes that have not yet been assigned
1210  for(unsigned n=0;n<n_node;n++)
1211  {
1212 
1213  //POSITIONAL EQUATAIONS
1214  //If the node is a hanging node
1215  if(node_pt(n)->is_hanging())
1216  {
1217  //Find the local hang info object
1218  HangInfo* hang_info_pt = node_pt(n)->hanging_pt();
1219  //Find the number of master nodes
1220  unsigned n_master = hang_info_pt->nmaster();
1221  //Loop over the master nodes
1222  for(unsigned m=0;m<n_master;m++)
1223  {
1224  //Get the m-th master node
1225  Node* Master_node_pt = hang_info_pt->master_node_pt(m);
1226 
1227  //If the local equation numbers associated with this master node
1228  //have not already been assigned, assign them
1229  if(local_eqn_number_done[Master_node_pt]==false)
1230  {
1231 
1232  //Now we need to test whether the master node is actually
1233  //a local node, in which case its local equation numbers
1234  //will already have been assigned and stored in
1235  //position_local_eqn(n,j,k)
1236 
1237  //Storage for the index of the local node
1238  //Initialised to n_node (beyond the possible indices of
1239  //"real" nodes)
1240  unsigned local_node_index=n_node;
1241  //Loop over the local nodes (again)
1242  for(unsigned n1=0;n1<n_node;n1++)
1243  {
1244  //If the master node is a local node
1245  //get its index and break out of the loop
1246  if(Master_node_pt==node_pt(n1))
1247  {
1248  local_node_index=n1;
1249  break;
1250  }
1251  }
1252 
1253  //Now we test whether the node was found
1254  if(local_node_index < n_node)
1255  {
1256  //Loop over the number of position dofs
1257  for(unsigned j=0;j<n_position_type;j++)
1258  {
1259  //Loop over the dimension of each node
1260  for(unsigned k=0;k<nodal_dim;k++)
1261  {
1262  //Set the values in the node-based positional look-up scheme
1263  Position_local_eqn_at_node(j,k) =
1264  position_local_eqn(local_node_index,j,k);
1265  }
1266  }
1267  }
1268  //Otherwise it's a new master node
1269  else
1270  {
1271  //Loop over the number of position dofs
1272  for(unsigned j=0;j<n_position_type;j++)
1273  {
1274  //Loop over the dimension of each node
1275  for(unsigned k=0;k<nodal_dim;k++)
1276  {
1277  //Get equation number (position_eqn_number)
1278  //Note eqn_number is long !
1279  long eqn_number =
1280  static_cast<SolidNode*>(Master_node_pt)
1281  ->position_eqn_number(j,k);
1282  //If equation_number positive add to array
1283  if(eqn_number >= 0)
1284  {
1285  //Add global equation number to the local queue
1286  global_eqn_number_queue.push_back(eqn_number);
1287  //Add pointer to the dof to the queue if required
1288  if(store_local_dof_pt)
1289  {
1291  &(Master_node_pt->x_gen(j,k)));
1292  }
1293  //Add to pointer-based scheme
1294  Position_local_eqn_at_node(j,k) = local_eqn_number;
1295  //Increase the number of local variables
1296  local_eqn_number++;
1297  }
1298  //Otherwise the value is pinned
1299  else
1300  {
1301  Position_local_eqn_at_node(j,k) = Data::Is_pinned;
1302  }
1303  }
1304  }
1305  } //End of case when it's a new master node
1306 
1307  //Dofs included with this node have now been done
1308  local_eqn_number_done[Master_node_pt] = true;
1309  //Add to the pointer-based reference scheme
1310  Local_position_hang_eqn[Master_node_pt] =
1311  Position_local_eqn_at_node;
1312  }
1313  }
1314  }
1315 
1316  } //End of loop over nodes
1317 
1318  //Now add our global equations numbers to the internal element storage
1319  add_global_eqn_numbers(global_eqn_number_queue,
1321  //Clear the memory used in the deque
1322  if(store_local_dof_pt)
1323  {std::deque<double*>().swap(GeneralisedElement::Dof_pt_deque);}
1324 
1325 
1326  } //End of if nodes
1327  }
1328 
1329 //============================================================================
1330 /// This function calculates the entries of Jacobian matrix, used in
1331 /// the Newton method, associated with the elastic problem in which the
1332 /// nodal position is a variable. It does this using finite differences,
1333 /// rather than an analytical formulation, so can be done in total generality.
1334 /// Overload the standard case to include hanging node case
1335 //==========================================================================
1338  DenseMatrix<double> &jacobian)
1339  {
1340  //Find the number of nodes
1341  const unsigned n_node = nnode();
1342 
1343  //If there are no nodes, return straight away
1344  if(n_node == 0) {return;}
1345 
1346  //Call the update function to ensure that the element is in
1347  //a consistent state before finite differencing starts
1349 
1350 // bool use_first_order_fd=false;
1351 
1352  //Find the number of positional dofs and nodal dimension
1353  const unsigned n_position_type = nnodal_position_type();
1354  const unsigned nodal_dim = nodal_dimension();
1355 
1356  //Find the number of dofs in the element
1357  const unsigned n_dof = ndof();
1358 
1359  //Create newres vector
1360  Vector<double> newres(n_dof); //, newres_minus(n_dof);
1361 
1362  //Used default value defined in GeneralisedElement
1363  const double fd_step = Default_fd_jacobian_step;
1364 
1365  //Integer storage for local unknowns
1366  int local_unknown=0;
1367 
1368  //Loop over the nodes
1369  for(unsigned l=0;l<n_node;l++)
1370  {
1371  //Get the pointer to the node
1372  Node* const local_node_pt = node_pt(l);
1373 
1374  //If the node is not a hanging node
1375  if(local_node_pt->is_hanging()==false)
1376  {
1377  //Loop over position dofs
1378  for(unsigned k=0;k<n_position_type;k++)
1379  {
1380  //Loop over dimension
1381  for(unsigned i=0;i<nodal_dim;i++)
1382  {
1383  local_unknown = position_local_eqn(l,k,i);
1384  //If the variable is free
1385  if(local_unknown >= 0)
1386  {
1387  //Store a pointer to the (generalised) Eulerian nodal position
1388  double* const value_pt = &(local_node_pt->x_gen(k,i));
1389 
1390  //Save the old value of the (generalised) Eulerian nodal position
1391  const double old_var = *value_pt;
1392 
1393  //Increment the (generalised) Eulerian nodal position
1394  *value_pt += fd_step;
1395 
1396  // Perform any auxialiary node updates
1397  local_node_pt->perform_auxiliary_node_update_fct();
1398 
1399  // Update any other slaved variables
1401 
1402 
1403  //Calculate the new residuals
1404  get_residuals(newres);
1405 
1406 // if (use_first_order_fd)
1407  {
1408  //Do forward finite differences
1409  for(unsigned m=0;m<n_dof;m++)
1410  {
1411  //Stick the entry into the Jacobian matrix
1412  jacobian(m,local_unknown) = (newres[m] - residuals[m])/fd_step;
1413  }
1414  }
1415 // else
1416 // {
1417 // //Take backwards step for the (generalised) Eulerian nodal
1418 // // position
1419 // node_pt(l)->x_gen(k,i) = old_var-fd_step;
1420 
1421 // //Calculate the new residuals at backward position
1422 // get_residuals(newres_minus);
1423 
1424 // //Do central finite differences
1425 // for(unsigned m=0;m<n_dof;m++)
1426 // {
1427 // //Stick the entry into the Jacobian matrix
1428 // jacobian(m,local_unknown) =
1429 // (newres[m] - newres_minus[m])/(2.0*fd_step);
1430 // }
1431 // }
1432 
1433  //Reset the (generalised) Eulerian nodal position
1434  *value_pt = old_var;
1435 
1436  // Perform any auxialiary node updates
1437  local_node_pt->perform_auxiliary_node_update_fct();
1438 
1439  // Reset any other slaved variables
1441  }
1442  }
1443  }
1444  }
1445  //Otherwise it's a hanging node
1446  else
1447  {
1448  //Find the local hanging object
1449  HangInfo* hang_info_pt = local_node_pt->hanging_pt();
1450  //Loop over the master nodes
1451  const unsigned n_master = hang_info_pt->nmaster();
1452  for(unsigned m=0;m<n_master;m++)
1453  {
1454  //Get the pointer to the master node
1455  Node* const master_node_pt = hang_info_pt->master_node_pt(m);
1456 
1457  //Get the local equation numbers for the master node
1458  DenseMatrix<int> Position_local_eqn_at_node
1459  = Local_position_hang_eqn[master_node_pt];
1460 
1461  //Loop over position dofs
1462  for(unsigned k=0;k<n_position_type;k++)
1463  {
1464  //Loop over dimension
1465  for(unsigned i=0;i<nodal_dim;i++)
1466  {
1467  local_unknown = Position_local_eqn_at_node(k,i);
1468  //If the variable is free
1469  if(local_unknown >= 0)
1470  {
1471  //Store a pointer to the (generalised) Eulerian nodal position
1472  double* const value_pt = &(master_node_pt->x_gen(k,i));
1473 
1474  //Save the old value of the (generalised) Eulerian nodal position
1475  const double old_var = *value_pt;
1476 
1477  //Increment the (generalised) Eulerian nodal position
1478  *value_pt += fd_step;
1479 
1480  // Perform any auxialiary node updates
1481  master_node_pt->perform_auxiliary_node_update_fct();
1482 
1483  // Update any slaved variables
1485 
1486  //Calculate the new residuals
1487  get_residuals(newres);
1488 
1489 // if (use_first_order_fd)
1490  {
1491  //Do forward finite differences
1492  for(unsigned m=0;m<n_dof;m++)
1493  {
1494  //Stick the entry into the Jacobian matrix
1495  jacobian(m,local_unknown) =
1496  (newres[m] - residuals[m])/fd_step;
1497  }
1498  }
1499 // else
1500 // {
1501 // //Take backwards step for the (generalised) Eulerian nodal
1502 // // position
1503 // master_node_pt->x_gen(k,i) = old_var-fd_step;
1504 
1505 // //Calculate the new residuals at backward position
1506 // get_residuals(newres_minus);
1507 
1508 // //Do central finite differences
1509 // for(unsigned m=0;m<n_dof;m++)
1510 // {
1511 // //Stick the entry into the Jacobian matrix
1512 // jacobian(m,local_unknown) =
1513 // (newres[m] - newres_minus[m])/(2.0*fd_step);
1514 // }
1515 // }
1516 
1517  //Reset the (generalised) Eulerian nodal position
1518  *value_pt = old_var;
1519 
1520  // Perform any auxialiary node updates
1521  master_node_pt->perform_auxiliary_node_update_fct();
1522 
1523  // Reset any other slaved variables
1525  }
1526  }
1527  }
1528  }
1529  } //End of hanging node case
1530 
1531  } //End of loop over nodes
1532 
1533  //End of finite difference loop
1534  //Final reset of any slaved data
1536  }
1537 
1538 }
int local_hang_eqn(Node *const &node_pt, const unsigned &i)
Access function that returns the local equation number for the hanging node variables (values stored ...
virtual void get_residuals(Vector< double > &residuals)
Calculate the vector of residuals of the equations in the element. By default initialise the vector t...
Definition: elements.h:975
virtual unsigned ncont_interpolated_values() const =0
Number of continuously interpolated values. Note: We assume that they are located at the beginning of...
virtual void reset_in_nodal_fd(const unsigned &i)
Function called within the finite difference loop for nodal data after the i-th nodal values is reset...
Definition: elements.h:1684
Data * geom_data_pt(const unsigned &j)
Return pointer to the j-th Data item that the object's shape depends on: Positional data of non-hangi...
void identify_field_data_for_interactions(std::set< std::pair< Data *, unsigned > > &paired_field_data)
The purpose of this function is to identify all possible.
void assemble_eulerian_base_vectors(const DShape &dpsids, DenseMatrix< double > &interpolated_G) const
Assemble the covariant Eulerian base vectors, assuming that the derivatives of the shape functions wi...
void assign_hanging_local_eqn_numbers(const bool &store_local_dof_pt)
Assign the local equation numbers for hanging node variables.
double nodal_position_gen(const unsigned &n, const unsigned &k, const unsigned &i) const
Return the value of the k-th type of the i-th positional variable at the local node n...
Definition: elements.h:2237
static double Max_integrity_tolerance
Max. allowed discrepancy in element integrity check.
cstr elem_len * i
Definition: cfortran.h:607
void get_dresidual_dnodal_coordinates(RankThreeTensor< double > &dresidual_dnodal_coordinates)
Compute derivatives of elemental residual vector with respect to nodal coordinates. Default implementation by FD can be overwritten for specific elements. dresidual_dnodal_coordinates(l,i,j) = d res(l) / dX_{ij} This version is overloaded from the version in FiniteElement and takes hanging nodes into account – j in the above loop loops over all the nodes that actively control the shape of the element (i.e. they are non-hanging or master nodes of hanging nodes in this element).
virtual void update_in_solid_position_fd(const unsigned &i)
Function called within the finite difference loop for the solid position dat after a change in any va...
Definition: elements.h:3997
static void check_value_id(const int &n_continuously_interpolated_values, const int &value_id)
Static helper function that is used to check that the value_id is in range.
int local_eqn_number(const unsigned long &ieqn_global) const
Return the local equation number corresponding to the ieqn_global-th global equation number...
Definition: elements.h:731
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
unsigned long eqn_number(const unsigned &ieqn_local) const
Return the global equation number corresponding to the ieqn_local-th local equation number...
Definition: elements.h:709
unsigned nnodal_position_type() const
Return the number of coordinate types that the element requires to interpolate the geometry between t...
Definition: elements.h:2340
double * value_pt(const unsigned &i) const
Return the pointer to the i-the stored value. Typically this is required when direct access to the st...
Definition: nodes.h:322
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
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 nodal_dimension() const
Return the required Eulerian dimension of the nodes in this element.
Definition: elements.h:2358
void add_global_eqn_numbers(std::deque< unsigned long > const &global_eqn_numbers, std::deque< double * > const &global_dof_pt)
Add the contents of the queue global_eqn_numbers to the local storage for the local-to-global transla...
Definition: elements.cc:149
long & eqn_number(const unsigned &i)
Return the equation number of the i-th stored variable.
Definition: nodes.h:365
void assemble_local_to_lagrangian_jacobian(const DShape &dpsids, DenseMatrix< double > &jacobian) const
Assemble the jacobian matrix for the mapping from local to lagrangian coordinates, given the derivatives of the shape function Overload the standard version to use the hanging information for the lagrangian coordinates.
unsigned nvalue() const
Return number of values stored in data object (incl pinned ones).
Definition: nodes.h:448
virtual void update_before_solid_position_fd()
Function that is called before the finite differencing of any solid position data. This may be overloaded to update any slaved data before finite differencing takes place.
Definition: elements.h:3987
void check_jacobian(const double &jacobian) const
Helper function used to check for singular or negative Jacobians in the transform from local to globa...
Definition: elements.cc:1722
unsigned nmaster() const
Return the number of master nodes.
Definition: nodes.h:733
std::map< Node *, int > * Local_hang_eqn
Storage for local equation numbers of hanging node variables (values stored at master nodes)...
unsigned ngeom_data() const
The number of geometric data affecting a RefineableSolidFiniteElement is the positional Data of all n...
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:995
bool is_hanging() const
Test whether the node is geometrically hanging.
Definition: nodes.h:1207
A Rank 3 Tensor class.
Definition: matrices.h:1337
virtual void update_before_nodal_fd()
Function that is called before the finite differencing of any nodal data. This may be overloaded to u...
Definition: elements.h:1670
static long Is_pinned
Static "Magic number" used in place of the equation number to indicate that the value is pinned...
Definition: nodes.h:192
virtual void fill_in_jacobian_from_nodal_by_fd(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Calculate the contributions to the jacobian from the nodal degrees of freedom using finite difference...
virtual void update_in_nodal_fd(const unsigned &i)
Function called within the finite difference loop for nodal data after a change in the i-th nodal val...
Definition: elements.h:1679
Data *& internal_data_pt(const unsigned &i)
Return a pointer to i-th internal data object.
Definition: elements.h:623
A class that represents a collection of data; each Data object may contain many different individual ...
Definition: nodes.h:89
double & x_gen(const unsigned &k, const unsigned &i)
Reference to the generalised position x(k,i).
Definition: nodes.h:1055
void assemble_local_to_eulerian_jacobian2(const DShape &d2psids, DenseMatrix< double > &jacobian2) const
Assemble the the "jacobian" matrix of second derivatives of the mapping from local to Eulerian coordi...
void fill_in_jacobian_from_solid_position_by_fd(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Compute element residual Vector and element Jacobian matrix corresponding to the solid positions...
static const unsigned N2deriv[]
Static array that holds the number of second derivatives as a function of the dimension of the elemen...
Definition: elements.h:1437
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
int nodal_local_eqn(const unsigned &n, const unsigned &i) const
Return the local equation number corresponding to the i-th value at the n-th local node...
Definition: elements.h:1383
Class that contains data for hanging nodes.
Definition: nodes.h:684
void assemble_local_to_eulerian_jacobian(const DShape &dpsids, DenseMatrix< double > &jacobian) const
Assemble the jacobian matrix for the mapping from local.
std::map< Node *, unsigned > Shape_controlling_node_lookup
Lookup scheme for unique number associated with any of the nodes that actively control the shape of t...
double lagrangian_position_gen(const unsigned &n, const unsigned &k, const unsigned &i) const
Return Generalised Lagrangian coordinate at local node n. `Direction' i, `Type' k.
Definition: elements.h:3646
unsigned ninternal_data() const
Return the number of internal data objects.
Definition: elements.h:828
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
unsigned nnodal_lagrangian_type() const
Return the number of types of (generalised) nodal Lagrangian coordinates required to interpolate the ...
Definition: elements.h:3526
unsigned ndof() const
Return the number of equations/dofs in the element.
Definition: elements.h:834
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition: elements.h:2470
double local_to_eulerian_mapping_diagonal(const DShape &dpsids, DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
Calculate the mapping from local to Eulerian coordinates given the derivatives of the shape functions...
std::map< Node *, DenseMatrix< int > > Local_position_hang_eqn
Storage for local equation numbers of hanging node variables associated with nodal positions...
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
virtual void reset_after_nodal_fd()
Function that is call after the finite differencing of the nodal data. This may be overloaded to rese...
Definition: elements.h:1675
static std::deque< double * > Dof_pt_deque
Static storage for deque used to add_global_equation_numbers when pointers to the dofs in each elemen...
Definition: elements.h:232
Node *const & master_node_pt(const unsigned &i) const
Return a pointer to the i-th master node.
Definition: nodes.h:736
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
Definition: nodes.h:992
virtual void deactivate_element()
Final operations that must be performed when the element is no longer active in the mesh...
void assign_solid_hanging_local_eqn_numbers(const bool &store_local_dof_pt)
Assign local equation numbers to the hanging values associated with positions or additional solid val...
virtual ~RefineableElement()
Destructor, delete the allocated storage for the hanging equations.
virtual void reset_in_solid_position_fd(const unsigned &i)
Function called within the finite difference loop for solid position data after the values in the i-t...
Definition: elements.h:4002
double local_to_lagrangian_mapping_diagonal(const DShape &dpsids, DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
Calculate the mapping from local to Lagrangian coordinates given the derivatives of the shape functio...
virtual void reset_after_solid_position_fd()
Function that is call after the finite differencing of the solid position data. This may be overloade...
Definition: elements.h:3992
int position_local_eqn(const unsigned &n, const unsigned &k, const unsigned &j) const
Access function that returns the local equation number that corresponds to the j-th coordinate of the...
Definition: elements.h:3881
static double Default_fd_jacobian_step
Double used for the default finite difference step in elemental jacobian calculations.
Definition: elements.h:1164
void assemble_local_to_lagrangian_jacobian2(const DShape &d2psids, DenseMatrix< double > &jacobian2) const
Assemble the the "jacobian" matrix of second derivatives, given the second derivatives of the shape f...
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...