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// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
7// LIC//
8// LIC// This library is free software; you can redistribute it and/or
9// LIC// modify it under the terms of the GNU Lesser General Public
10// LIC// License as published by the Free Software Foundation; either
11// LIC// version 2.1 of the License, or (at your option) any later version.
12// LIC//
13// LIC// This library is distributed in the hope that it will be useful,
14// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16// LIC// Lesser General Public License for more details.
17// LIC//
18// LIC// You should have received a copy of the GNU Lesser General Public
19// LIC// License along with this library; if not, write to the Free Software
20// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21// LIC// 02110-1301 USA.
22// LIC//
23// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24// LIC//
25// LIC//====================================================================
26// Non-inline member functions for refineable elements
27
28// oomph-lib includes
29#include "refineable_elements.h"
30#include "shape.h"
31
32namespace oomph
33{
34 //=====================================================================
35 /// Destructor (needed here because of IBM xlC compiler under AIX)
36 /// Delete the storage allocated for the local equations.
37 //====================================================================
39 {
40 delete[] Local_hang_eqn;
41 }
42
43 //========================================================================
44 /// Max. allowed discrepancy in element integrity check
45 //========================================================================
47
48 //=========================================================================
49 /// Helper function that is used to check that the value_id is in the range
50 /// allowed by the element. The number of continuously interpolated values
51 /// and the value_id must be passed as arguments.
52 //========================================================================
54 const int& n_continuously_interpolated_values, const int& value_id)
55 {
56 // If the value_id is more than the number of continuously interpolated
57 // values or less than -1 (for the position), we have a problem
58 if ((value_id < -1) || (value_id >= n_continuously_interpolated_values))
59 {
60 std::ostringstream error_stream;
61 error_stream << "Value_id " << value_id << " is out of range."
62 << std::endl
63 << "It can only take the values -1 (position) "
64 << "or an integer in the range 0 to "
65 << n_continuously_interpolated_values << std::endl;
66
67 throw OomphLibError(
68 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
69 }
70 }
71
72 //=========================================================================
73 /// Internal function that is used to assemble the jacobian of the mapping
74 /// from local coordinates (s) to the eulerian coordinates (x), given the
75 /// derivatives of the shape functions.
76 //=========================================================================
78 const DShape& dpsids, DenseMatrix<double>& jacobian) const
79 {
80 // Locally cache the elemental dimension
81 const unsigned el_dim = dim();
82 // The number of shape functions must be equal to the number
83 // of nodes (by definition)
84 const unsigned n_shape = nnode();
85 // The number of shape function types must be equal to the number
86 // of nodal position types (by definition)
87 const unsigned n_shape_type = nnodal_position_type();
88
89#ifdef PARANOID
90 // Check for dimensional compatibility
91 if (dim() != nodal_dimension())
92 {
93 std::ostringstream error_message;
94 error_message << "Dimension mismatch" << std::endl;
95 error_message << "The elemental dimension: " << dim()
96 << " must equal the nodal dimension: " << nodal_dimension()
97 << " for the jacobian of the mapping to be well-defined"
98 << std::endl;
99 throw OomphLibError(
100 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
101 }
102#endif
103
104 // Loop over the rows of the jacobian
105 for (unsigned i = 0; i < el_dim; i++)
106 {
107 // Loop over the columns of the jacobian
108 for (unsigned j = 0; j < el_dim; j++)
109 {
110 // Zero the entry
111 jacobian(i, j) = 0.0;
112 // Loop over the shape functions
113 for (unsigned l = 0; l < n_shape; l++)
114 {
115 for (unsigned k = 0; k < n_shape_type; k++)
116 {
117 // Jacobian is dx_j/ds_i, which is represented by the sum
118 // over the dpsi/ds_i of the nodal points X j
119 // Call the Hanging version of positions
120 jacobian(i, j) += nodal_position_gen(l, k, j) * dpsids(l, k, i);
121 }
122 }
123 }
124 }
125 }
126
127 //=========================================================================
128 /// Internal function that is used to assemble the jacobian of second
129 /// derivatives of the the mapping from local coordinates (s) to the
130 /// eulerian coordinates (x), given the second derivatives of the
131 /// shape functions.
132 //=========================================================================
134 const DShape& d2psids, DenseMatrix<double>& jacobian2) const
135 {
136 // Find the the dimension of the element
137 const unsigned el_dim = dim();
138 // Find the number of shape functions and shape functions types
139 // Must be equal to the number of nodes and their position types
140 // by the definition of the shape function.
141 const unsigned n_shape = nnode();
142 const unsigned n_shape_type = nnodal_position_type();
143 // Find the number of second derivatives
144 const unsigned n_row = N2deriv[el_dim];
145
146 // Assemble the "jacobian" (d^2 x_j/ds_i^2) for second derivatives of
147 // shape functions
148 // Loop over the rows (number of second derivatives)
149 for (unsigned i = 0; i < n_row; i++)
150 {
151 // Loop over the columns (element dimension
152 for (unsigned j = 0; j < el_dim; j++)
153 {
154 // Zero the entry
155 jacobian2(i, j) = 0.0;
156 // Loop over the shape functions
157 for (unsigned l = 0; l < n_shape; l++)
158 {
159 // Loop over the shape function types
160 for (unsigned k = 0; k < n_shape_type; k++)
161 {
162 // Add the terms to the jacobian entry
163 // Call the Hanging version of positions
164 jacobian2(i, j) += nodal_position_gen(l, k, j) * d2psids(l, k, i);
165 }
166 }
167 }
168 }
169 }
170
171 //=====================================================================
172 /// Assemble the covariant Eulerian base vectors and return them in the
173 /// matrix interpolated_G. The derivatives of the shape functions with
174 /// respect to the local coordinate should already have been calculated
175 /// before calling this function
176 //=====================================================================
178 const DShape& dpsids, DenseMatrix<double>& interpolated_G) const
179 {
180 // Find the number of nodes and position types
181 const unsigned n_node = nnode();
182 const unsigned n_position_type = nnodal_position_type();
183 // Find the dimension of the node and element
184 const unsigned n_dim_node = nodal_dimension();
185 const unsigned n_dim_element = dim();
186
187 // Loop over the dimensions and compute the entries of the
188 // base vector matrix
189 for (unsigned i = 0; i < n_dim_element; i++)
190 {
191 for (unsigned j = 0; j < n_dim_node; j++)
192 {
193 // Initialise the j-th component of the i-th base vector to zero
194 interpolated_G(i, j) = 0.0;
195 for (unsigned l = 0; l < n_node; l++)
196 {
197 for (unsigned k = 0; k < n_position_type; k++)
198 {
199 interpolated_G(i, j) +=
200 nodal_position_gen(l, k, j) * dpsids(l, k, i);
201 }
202 }
203 }
204 }
205 }
206
207
208 //==========================================================================
209 /// Calculate the mapping from local to eulerian coordinates
210 /// assuming that the coordinates are aligned in the direction of the local
211 /// coordinates, i.e. there are no cross terms and the jacobian is diagonal.
212 /// The local derivatives are passed as dpsids and the jacobian and
213 /// inverse jacobian are returned.
214 //==========================================================================
216 const DShape& dpsids,
217 DenseMatrix<double>& jacobian,
218 DenseMatrix<double>& inverse_jacobian) const
219 {
220 // Find the dimension of the element
221 const unsigned el_dim = dim();
222 // Find the number of shape functions and shape functions types
223 // Equal to the number of nodes and their position types by definition
224 const unsigned n_shape = nnode();
225 const unsigned n_shape_type = nnodal_position_type();
226
227#ifdef PARANOID
228 // Check for dimension compatibility
229 if (dim() != nodal_dimension())
230 {
231 std::ostringstream error_message;
232 error_message << "Dimension mismatch" << std::endl;
233 error_message << "The elemental dimension: " << dim()
234 << " must equal the nodal dimension: " << nodal_dimension()
235 << " for the jacobian of the mapping to be well-defined"
236 << std::endl;
237 throw OomphLibError(
238 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
239 }
240#endif
241
242 // In this case we assume that there are no cross terms, that is
243 // global coordinate 0 is always in the direction of local coordinate 0
244
245 // Loop over the coordinates
246 for (unsigned i = 0; i < el_dim; i++)
247 {
248 // Zero the jacobian and inverse jacobian entries
249 for (unsigned j = 0; j < el_dim; j++)
250 {
251 jacobian(i, j) = 0.0;
252 inverse_jacobian(i, j) = 0.0;
253 }
254
255 // Loop over the shape functions
256 for (unsigned l = 0; l < n_shape; l++)
257 {
258 // Loop over the types of dof
259 for (unsigned k = 0; k < n_shape_type; k++)
260 {
261 // Derivatives are always dx_{i}/ds_{i}
262 jacobian(i, i) += nodal_position_gen(l, k, i) * dpsids(l, k, i);
263 }
264 }
265 }
266
267 // Now calculate the determinant of the matrix
268 double det = 1.0;
269 for (unsigned i = 0; i < el_dim; i++)
270 {
271 det *= jacobian(i, i);
272 }
273
274// Report if Matrix is singular, or negative
275#ifdef PARANOID
276 check_jacobian(det);
277#endif
278
279 // Calculate the inverse mapping (trivial in this case)
280 for (unsigned i = 0; i < el_dim; i++)
281 {
282 inverse_jacobian(i, i) = 1.0 / jacobian(i, i);
283 }
284
285 // Return the value of the Jacobian
286 return (det);
287 }
288
289 //========================================================================
290 /// Deactivate the element by marking setting all local pointers to
291 /// obsolete nodes to zero
292 //=======================================================================
294 {
295 // Find the number of nodes
296 const unsigned n_node = nnode();
297 // Loop over the nodes
298 for (unsigned n = 0; n < n_node; n++)
299 {
300 // If the node pointer has not already been nulled, but is obsolete
301 // then null it
302 if ((this->node_pt(n) != 0) && (this->node_pt(n)->is_obsolete()))
303 {
304 this->node_pt(n) = 0;
305 }
306 }
307 }
308
309 //=======================================================================
310 /// Assign the local hang eqn.
311 //=======================================================================
313 const bool& store_local_dof_pt)
314 {
315 // Find the number of nodes
316 const unsigned n_node = nnode();
317
318 // Check there are nodes!
319 if (n_node > 0)
320 {
321 // Find the number of continuously interpolated values
322 const unsigned n_cont_values = ncont_interpolated_values();
323
324 // Delete the storage associated with the local hanging equation schemes
325 delete[] Local_hang_eqn;
326 // Allocate new storage
327 Local_hang_eqn = new std::map<Node*, int>[n_cont_values];
328
329 // Boolean that is set to true if there are hanging equation numbers
330 bool hanging_eqn_numbers = false;
331
332 // Maps that store whether the node's local equation number for a
333 // particular value has already been assigned
334 Vector<std::map<Node*, bool>> local_eqn_number_done(n_cont_values);
335
336 // Get number of dofs assigned thus far
337 unsigned local_eqn_number = ndof();
338
339 // A local queue to store the global equation numbers
340 std::deque<unsigned long> global_eqn_number_queue;
341
342 // Now loop over all the nodes again to find the master nodes
343 // external to the element
344 for (unsigned n = 0; n < n_node; n++)
345 {
346 // Loop over the number of continuously-interpolated values at the node
347 for (unsigned j = 0; j < n_cont_values; j++)
348 {
349 // If the node is hanging node in value j
350 if (node_pt(n)->is_hanging(j))
351 {
352 // Get the pointer to the appropriate hanging info object
353 HangInfo* hang_info_pt = node_pt(n)->hanging_pt(j);
354 // Find the number of master nodes
355 unsigned n_master = hang_info_pt->nmaster();
356 // Loop over the master nodes
357 for (unsigned m = 0; m < n_master; m++)
358 {
359 // Get the m-th master node
360 Node* Master_node_pt = hang_info_pt->master_node_pt(m);
361
362 // If the master node's value has not been considered already,
363 // give it a local equation number
364 if (local_eqn_number_done[j][Master_node_pt] == false)
365 {
366#ifdef PARANOID
367 // Check that the value is stored at the master node
368 // If not the master node must have be set up incorrectly
369 if (Master_node_pt->nvalue() < j)
370 {
371 std::ostringstream error_stream;
372 error_stream << "Master node for " << j
373 << "-th value only has "
374 << Master_node_pt->nvalue() << " stored values!"
375 << std::endl;
376
377 throw OomphLibError(error_stream.str(),
378 OOMPH_CURRENT_FUNCTION,
379 OOMPH_EXCEPTION_LOCATION);
380 }
381#endif
382
383 // We now need to test whether the master node is actually
384 // a local node, in which case its local equation number
385 // will already have been assigned and will be stored
386 // in nodal_local_eqn
387
388 // Storage for the index of the local node
389 // Initialised to n_node (beyond the possible indices of
390 //"real" nodes)
391 unsigned local_node_index = n_node;
392 // Loop over the local nodes (again)
393 for (unsigned n1 = 0; n1 < n_node; n1++)
394 {
395 // If the master node is a local node
396 // get its index and break out of the loop
397 if (Master_node_pt == node_pt(n1))
398 {
399 local_node_index = n1;
400 break;
401 }
402 }
403
404 // Now we test whether the node was found
405 if (local_node_index < n_node)
406 {
407 // Copy the local equation number to the
408 // pointer-based look-up scheme
409 Local_hang_eqn[j][Master_node_pt] =
410 nodal_local_eqn(local_node_index, j);
411 }
412 // Otherwise it's a new master node
413 else
414 {
415 // Get the equation number
416 long eqn_number = Master_node_pt->eqn_number(j);
417 // If equation_number positive add to array
418 if (eqn_number >= 0)
419 {
420 // Add global equation number to the queue
421 global_eqn_number_queue.push_back(eqn_number);
422 // Add pointer to the dof to the queue if required
423 if (store_local_dof_pt)
424 {
426 Master_node_pt->value_pt(j));
427 }
428 // Add to pointer based scheme
429 Local_hang_eqn[j][Master_node_pt] = local_eqn_number;
430 // Increase number of local variables
432 }
433 // Otherwise the value is pinned
434 else
435 {
436 Local_hang_eqn[j][Master_node_pt] = Data::Is_pinned;
437 }
438 // There are now hanging equation numbers
439 }
440 // The value at this node has now been done
441 local_eqn_number_done[j][Master_node_pt] = true;
442 // There are hanging equation numbers
443 hanging_eqn_numbers = true;
444 }
445 }
446 }
447 }
448 } // End of second loop over nodes
449
450 // Now add our global equations numbers to the internal element storage
451 add_global_eqn_numbers(global_eqn_number_queue,
453 // Clear the memory used in the deque
454 if (store_local_dof_pt)
455 {
456 std::deque<double*>().swap(GeneralisedElement::Dof_pt_deque);
457 }
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 /// Compute derivatives of elemental residual vector with respect
588 /// to nodal coordinates. Default implementation by FD can be overwritten
589 /// for specific elements.
590 /// dresidual_dnodal_coordinates(l,i,j) = d res(l) / dX_{ij}
591 /// This version is overloaded from the version in FiniteElement
592 /// and takes hanging nodes into account -- j in the above loop
593 /// loops over all the nodes that actively control the
594 /// shape of the element (i.e. they are non-hanging or master nodes of
595 /// hanging nodes in this element).
596 //==========================================================================
598 RankThreeTensor<double>& dresidual_dnodal_coordinates)
599 {
600 // Number of nodes
601 unsigned n_nod = nnode();
602
603 // If the element has no nodes (why??!!) return straightaway
604 if (n_nod == 0) return;
605
606 // Get dimension from first node
607 unsigned dim_nod = node_pt(0)->ndim();
608
609 // Number of dofs
610 unsigned n_dof = ndof();
611
612 // Get reference residual
613 Vector<double> res(n_dof);
614 Vector<double> res_pls(n_dof);
615 get_residuals(res);
616
617 // FD step
619
620 // Do FD loop over all active nodes
621 for (std::map<Node*, unsigned>::iterator it =
624 it++)
625 {
626 // Get node
627 Node* nod_pt = it->first;
628
629 // Get its number
630 unsigned node_number = it->second;
631
632 // Loop over coordinate directions
633 for (unsigned i = 0; i < dim_nod; i++)
634 {
635 // Make backup
636 double backup = nod_pt->x(i);
637
638 // Do FD step. No node update required as we're
639 // attacking the coordinate directly...
640 nod_pt->x(i) += eps_fd;
641
642 // Perform auxiliary node update function
644
645 // Get advanced residual
646 get_residuals(res_pls);
647
648 // Fill in FD entries [Loop order is "wrong" here as l is the
649 // slow index but this is in a function that's costly anyway
650 // and gives us the fastest loop outside where these tensor
651 // is actually used.]
652 for (unsigned l = 0; l < n_dof; l++)
653 {
654 dresidual_dnodal_coordinates(l, i, node_number) =
655 (res_pls[l] - res[l]) / eps_fd;
656 }
657
658 // Reset coordinate. No node update required as we're
659 // attacking the coordinate directly...
660 nod_pt->x(i) = backup;
661
662 // Perform auxiliary node update function
664 }
665 }
666 }
667
668
669 //============================================================================
670 /// This function calculates the entries of Jacobian matrix, used in
671 /// the Newton method, associated with the nodal degrees of freedom.
672 /// This is done by finite differences to handle the general case
673 /// Overload the standard case to include hanging node case
674 //==========================================================================
676 Vector<double>& residuals, DenseMatrix<double>& jacobian)
677 {
678 // Find the number of nodes
679 const unsigned n_node = nnode();
680
681 // If there are no nodes, return straight away
682 if (n_node == 0)
683 {
684 return;
685 }
686
687 // Call the update function to ensure that the element is in
688 // a consistent state before finite differencing starts
690
691 // bool use_first_order_fd=false;
692
693 // Find the number of dofs in the element
694 const unsigned n_dof = ndof();
695
696 // Create newres vector
697 Vector<double> newres(n_dof); // , newres_minus(n_dof);
698
699 // Used default value defined in GeneralisedElement
700 const double fd_step = Default_fd_jacobian_step;
701
702 // Integer storage for local unknowns
703 int local_unknown = 0;
704
705 // Loop over the nodes
706 for (unsigned n = 0; n < n_node; n++)
707 {
708 // Get the number of values stored at the node
709 const unsigned n_value = node_pt(n)->nvalue();
710
711 // Get a pointer to the local node
712 Node* const local_node_pt = node_pt(n);
713
714 // Loop over the number of values stored at the node
715 for (unsigned i = 0; i < n_value; i++)
716 {
717 // If the nodal value is not hanging
718 if (local_node_pt->is_hanging(i) == false)
719 {
720 // Get the equation number
721 local_unknown = nodal_local_eqn(n, i);
722 // If the variable is free
723 if (local_unknown >= 0)
724 {
725 // Get a pointer to the nodal value
726 double* const value_pt = local_node_pt->value_pt(i);
727
728 // Save the old value of nodal value
729 const double old_var = *value_pt;
730
731 // Increment the nodal value
732 *value_pt += fd_step;
733
734 // Now update any dependent variables
736
737 // Calculate the new residuals
738 get_residuals(newres);
739
740 // if (use_first_order_fd)
741 {
742 // Do forward finite differences
743 for (unsigned m = 0; m < n_dof; m++)
744 {
745 // Stick the entry into the Jacobian matrix
746 jacobian(m, local_unknown) =
747 (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 dependent 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 dependent 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) =
813 (newres[m] - residuals[m]) / fd_step;
814 }
815 }
816 // else
817 // {
818 // //Take backwards step
819 // *value_pt = old_var - fd_step;
820
821 // //Calculate the new residuals at backward
822 // position get_residuals(newres_minus);
823
824 // //Do central finite differences
825 // for(unsigned m=0;m<n_dof;m++)
826 // {
827 // //Stick the entry into the Jacobian matrix
828 // jacobian(m,local_unknown) =
829 // (newres[m] - newres_minus[m])/(2.0*fd_step);
830 // }
831 // }
832
833 // Reset the value at the master node
834 *value_pt = old_var;
835
836 // Reset any dependent variables
838 }
839 } // End of loop over master nodes
840 } // End of hanging node case
841 } // End of loop over values
842 } // End of loop over nodes
843
844 // End of finite difference loop
845 // Final reset of any dependent data
847 }
848
849
850 //=========================================================================
851 /// Internal function that is used to assemble the jacobian of the mapping
852 /// from local coordinates (s) to the lagrangian coordinates (xi), given the
853 /// derivatives of the shape functions.
854 //=========================================================================
856 const DShape& dpsids, 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) +=
881 lagrangian_position_gen(l, k, j) * dpsids(l, k, i);
882 }
883 }
884 }
885 }
886 }
887
888 //=========================================================================
889 /// Internal function that is used to assemble the jacobian of second
890 /// derivatives of the the mapping from local coordinates (s) to the
891 /// lagrangian coordinates (xi), given the second derivatives of the
892 /// shape functions.
893 //=========================================================================
895 const DShape& d2psids, 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) +=
924 lagrangian_position_gen(l, k, j) * d2psids(l, k, i);
925 }
926 }
927 }
928 }
929 }
930
931 //==========================================================================
932 /// Calculate the mapping from local to lagrangian coordinates
933 /// assuming that the coordinates are aligned in the direction of the local
934 /// coordinates, i.e. there are no cross terms and the jacobian is diagonal.
935 /// The local derivatives are passed as dpsids and the jacobian and
936 /// inverse jacobian are returned.
937 //==========================================================================
939 const DShape& dpsids,
940 DenseMatrix<double>& jacobian,
941 DenseMatrix<double>& inverse_jacobian) 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 {
958 jacobian(i, j) = 0.0;
959 inverse_jacobian(i, j) = 0.0;
960 }
961
962 // Loop over the shape functions
963 for (unsigned l = 0; l < n_shape; l++)
964 {
965 // Loop over the types of dof
966 for (unsigned k = 0; k < n_shape_type; k++)
967 {
968 // Derivatives are always dx_{i}/ds_{i}
969 jacobian(i, i) += lagrangian_position_gen(l, k, i) * dpsids(l, k, i);
970 }
971 }
972 }
973
974 // Now calculate the determinant of the matrix
975 double det = 1.0;
976 for (unsigned i = 0; i < el_dim; i++)
977 {
978 det *= jacobian(i, i);
979 }
980
981// Report if Matrix is singular, or negative
982#ifdef PARANOID
983 check_jacobian(det);
984#endif
985
986 // Calculate the inverse mapping (trivial in this case)
987 for (unsigned i = 0; i < el_dim; i++)
988 {
989 inverse_jacobian(i, i) = 1.0 / jacobian(i, i);
990 }
991
992 // Return the value of the Jacobian
993 return (det);
994 }
995
996
997 //========================================================================
998 /// The number of geometric data affecting a
999 /// RefineableSolidFiniteElement is the positional Data of all
1000 /// non-hanging nodes plus the geometric Data of all distinct
1001 /// master nodes. Recomputed on the fly.
1002 //========================================================================
1004 {
1005 // Find the number of nodes
1006 const unsigned n_node = nnode();
1007
1008 // Temporary storage for unique position data
1009 std::set<Data*> all_position_data_pt;
1010
1011 // Now loop over all the nodes again to find the master nodes
1012 // of any hanging nodes that have not yet been assigned
1013 for (unsigned n = 0; n < n_node; n++)
1014 {
1015 // If the node is a hanging node
1016 if (node_pt(n)->is_hanging())
1017 {
1018 // Find the local hang info object
1019 HangInfo* hang_info_pt = node_pt(n)->hanging_pt();
1020
1021 // Find the number of master nodes
1022 unsigned n_master = hang_info_pt->nmaster();
1023
1024 // Loop over the master nodes
1025 for (unsigned m = 0; m < n_master; m++)
1026 {
1027 // Get the m-th master node
1028 Node* Master_node_pt = hang_info_pt->master_node_pt(m);
1029
1030 // Add to set
1031 all_position_data_pt.insert(
1032 dynamic_cast<SolidNode*>(Master_node_pt)->variable_position_pt());
1033 }
1034 }
1035 // Not hanging
1036 else
1037 {
1038 // Add node itself to set
1039 all_position_data_pt.insert(
1040 dynamic_cast<SolidNode*>(node_pt(n))->variable_position_pt());
1041 }
1042
1043 } // End of loop over nodes
1044
1045 // How many are there?
1046 return all_position_data_pt.size();
1047 }
1048
1049
1050 //========================================================================
1051 /// Return pointer to the j-th Data item that the object's
1052 /// shape depends on: Positional data of non-hanging nodes and
1053 /// positional data of master nodes. Recomputed on the fly.
1054 //========================================================================
1056 {
1057 // Find the number of nodes
1058 const unsigned n_node = nnode();
1059
1060 // Temporary storage for unique position data. Set and vector are
1061 // required to ensure uniqueness in enumeration on different processors.
1062 // Set checks uniqueness; Vector stores entries in predictable order.
1063 std::set<Data*> all_position_data_pt;
1064 Vector<Data*> all_position_data_vector_pt;
1065
1066 // Number of entries in set before possibly adding new entry
1067 unsigned n_old = 0;
1068
1069 // Now loop over all the nodes again to find the master nodes
1070 // of any hanging nodes that have not yet been assigned
1071 for (unsigned n = 0; n < n_node; n++)
1072 {
1073 // If the node is a hanging node
1074 if (node_pt(n)->is_hanging())
1075 {
1076 // Find the local hang info object
1077 HangInfo* hang_info_pt = node_pt(n)->hanging_pt();
1078
1079 // Find the number of master nodes
1080 unsigned n_master = hang_info_pt->nmaster();
1081
1082 // Loop over the master nodes
1083 for (unsigned m = 0; m < n_master; m++)
1084 {
1085 // Get the m-th master node
1086 Node* Master_node_pt = hang_info_pt->master_node_pt(m);
1087
1088 // Positional data
1089 Data* pos_data_pt =
1090 dynamic_cast<SolidNode*>(Master_node_pt)->variable_position_pt();
1091
1092 // Add to set
1093 n_old = all_position_data_pt.size();
1094 all_position_data_pt.insert(pos_data_pt);
1095
1096 // New entry?
1097 if (all_position_data_pt.size() > n_old)
1098 {
1099 all_position_data_vector_pt.push_back(pos_data_pt);
1100 }
1101 }
1102 }
1103 // Not hanging
1104 else
1105 {
1106 // Add node itself to set
1107
1108 // Positional data
1109 Data* pos_data_pt =
1110 dynamic_cast<SolidNode*>(node_pt(n))->variable_position_pt();
1111
1112 // Add to set
1113 n_old = all_position_data_pt.size();
1114 all_position_data_pt.insert(pos_data_pt);
1115
1116 // New entry?
1117 if (all_position_data_pt.size() > n_old)
1118 {
1119 all_position_data_vector_pt.push_back(pos_data_pt);
1120 }
1121 }
1122
1123 } // End of loop over nodes
1124
1125
1126 // Return j-th entry
1127 return all_position_data_vector_pt[j];
1128 }
1129
1130
1131 //========================================================================
1132 /// Specify Data that affects the geometry of the element
1133 /// by adding the position Data to the set that's passed in.
1134 /// (This functionality is required in FSI problems; set is used to
1135 /// avoid double counting). Refineable version includes hanging nodes
1136 //========================================================================
1138 std::set<Data*>& geometric_data_pt)
1139 {
1140 // Loop over the node update data and add to the set
1141 const unsigned n_node = this->nnode();
1142 for (unsigned j = 0; j < n_node; j++)
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
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, nodal_dim);
1196
1197 // Map that store whether the node's equation numbers have already been
1198 // added to the local arrays
1199 std::map<Node*, bool> local_eqn_number_done;
1200
1201 // Get number of dofs so far
1202 unsigned local_eqn_number = ndof();
1203
1204 // A local queue to store the global equation numbers
1205 std::deque<unsigned long> global_eqn_number_queue;
1206
1207 // Now loop over all the nodes again to find the master nodes
1208 // of any hanging nodes that have not yet been assigned
1209 for (unsigned n = 0; n < n_node; n++)
1210 {
1211 // POSITIONAL EQUATAIONS
1212 // If the node is a hanging node
1213 if (node_pt(n)->is_hanging())
1214 {
1215 // Find the local hang info object
1216 HangInfo* hang_info_pt = node_pt(n)->hanging_pt();
1217 // Find the number of master nodes
1218 unsigned n_master = hang_info_pt->nmaster();
1219 // Loop over the master nodes
1220 for (unsigned m = 0; m < n_master; m++)
1221 {
1222 // Get the m-th master node
1223 Node* Master_node_pt = hang_info_pt->master_node_pt(m);
1224
1225 // If the local equation numbers associated with this master node
1226 // have not already been assigned, assign them
1227 if (local_eqn_number_done[Master_node_pt] == false)
1228 {
1229 // Now we need to test whether the master node is actually
1230 // a local node, in which case its local equation numbers
1231 // will already have been assigned and stored in
1232 // position_local_eqn(n,j,k)
1233
1234 // Storage for the index of the local node
1235 // Initialised to n_node (beyond the possible indices of
1236 //"real" nodes)
1237 unsigned local_node_index = n_node;
1238 // Loop over the local nodes (again)
1239 for (unsigned n1 = 0; n1 < n_node; n1++)
1240 {
1241 // If the master node is a local node
1242 // get its index and break out of the loop
1243 if (Master_node_pt == node_pt(n1))
1244 {
1245 local_node_index = n1;
1246 break;
1247 }
1248 }
1249
1250 // Now we test whether the node was found
1251 if (local_node_index < n_node)
1252 {
1253 // Loop over the number of position dofs
1254 for (unsigned j = 0; j < n_position_type; j++)
1255 {
1256 // Loop over the dimension of each node
1257 for (unsigned k = 0; k < nodal_dim; k++)
1258 {
1259 // Set the values in the node-based positional look-up
1260 // scheme
1261 Position_local_eqn_at_node(j, k) =
1262 position_local_eqn(local_node_index, j, k);
1263 }
1264 }
1265 }
1266 // Otherwise it's a new master node
1267 else
1268 {
1269 // Loop over the number of position dofs
1270 for (unsigned j = 0; j < n_position_type; j++)
1271 {
1272 // Loop over the dimension of each node
1273 for (unsigned k = 0; k < nodal_dim; k++)
1274 {
1275 // Get equation number (position_eqn_number)
1276 // Note eqn_number is long !
1277 long eqn_number = static_cast<SolidNode*>(Master_node_pt)
1278 ->position_eqn_number(j, k);
1279 // If equation_number positive add to array
1280 if (eqn_number >= 0)
1281 {
1282 // Add global equation number to the local queue
1283 global_eqn_number_queue.push_back(eqn_number);
1284 // Add pointer to the dof to the queue if required
1285 if (store_local_dof_pt)
1286 {
1288 &(Master_node_pt->x_gen(j, k)));
1289 }
1290 // Add to pointer-based scheme
1291 Position_local_eqn_at_node(j, k) = local_eqn_number;
1292 // Increase the number of local variables
1294 }
1295 // Otherwise the value is pinned
1296 else
1297 {
1298 Position_local_eqn_at_node(j, k) = Data::Is_pinned;
1299 }
1300 }
1301 }
1302 } // End of case when it's a new master node
1303
1304 // Dofs included with this node have now been done
1305 local_eqn_number_done[Master_node_pt] = true;
1306 // Add to the pointer-based reference scheme
1307 Local_position_hang_eqn[Master_node_pt] =
1308 Position_local_eqn_at_node;
1309 }
1310 }
1311 }
1312
1313 } // End of loop over nodes
1314
1315 // Now add our global equations numbers to the internal element storage
1316 add_global_eqn_numbers(global_eqn_number_queue,
1318 // Clear the memory used in the deque
1319 if (store_local_dof_pt)
1320 {
1321 std::deque<double*>().swap(GeneralisedElement::Dof_pt_deque);
1322 }
1323
1324
1325 } // End of if nodes
1326 }
1327
1328 //============================================================================
1329 /// This function calculates the entries of Jacobian matrix, used in
1330 /// the Newton method, associated with the elastic problem in which the
1331 /// nodal position is a variable. It does this using finite differences,
1332 /// rather than an analytical formulation, so can be done in total generality.
1333 /// Overload the standard case to include hanging node case
1334 //==========================================================================
1336 Vector<double>& residuals, DenseMatrix<double>& jacobian)
1337 {
1338 // Find the number of nodes
1339 const unsigned n_node = nnode();
1340
1341 // If there are no nodes, return straight away
1342 if (n_node == 0)
1343 {
1344 return;
1345 }
1346
1347 // Call the update function to ensure that the element is in
1348 // a consistent state before finite differencing starts
1350
1351 // bool use_first_order_fd=false;
1352
1353 // Find the number of positional dofs and nodal dimension
1354 const unsigned n_position_type = nnodal_position_type();
1355 const unsigned nodal_dim = nodal_dimension();
1356
1357 // Find the number of dofs in the element
1358 const unsigned n_dof = ndof();
1359
1360 // Create newres vector
1361 Vector<double> newres(n_dof); //, newres_minus(n_dof);
1362
1363 // Used default value defined in GeneralisedElement
1364 const double fd_step = Default_fd_jacobian_step;
1365
1366 // Integer storage for local unknowns
1367 int local_unknown = 0;
1368
1369 // Loop over the nodes
1370 for (unsigned l = 0; l < n_node; l++)
1371 {
1372 // Get the pointer to the node
1373 Node* const local_node_pt = node_pt(l);
1374
1375 // If the node is not a hanging node
1376 if (local_node_pt->is_hanging() == false)
1377 {
1378 // Loop over position dofs
1379 for (unsigned k = 0; k < n_position_type; k++)
1380 {
1381 // Loop over dimension
1382 for (unsigned i = 0; i < nodal_dim; i++)
1383 {
1384 local_unknown = position_local_eqn(l, k, i);
1385 // If the variable is free
1386 if (local_unknown >= 0)
1387 {
1388 // Store a pointer to the (generalised) Eulerian nodal position
1389 double* const value_pt = &(local_node_pt->x_gen(k, i));
1390
1391 // Save the old value of the (generalised) Eulerian nodal position
1392 const double old_var = *value_pt;
1393
1394 // Increment the (generalised) Eulerian nodal position
1395 *value_pt += fd_step;
1396
1397 // Perform any auxialiary node updates
1398 local_node_pt->perform_auxiliary_node_update_fct();
1399
1400 // Update any other dependent variables
1402
1403
1404 // Calculate the new residuals
1405 get_residuals(newres);
1406
1407 // if (use_first_order_fd)
1408 {
1409 // Do forward finite differences
1410 for (unsigned m = 0; m < n_dof; m++)
1411 {
1412 // Stick the entry into the Jacobian matrix
1413 jacobian(m, local_unknown) =
1414 (newres[m] - residuals[m]) / fd_step;
1415 }
1416 }
1417 // else
1418 // {
1419 // //Take backwards step for the (generalised)
1420 // Eulerian nodal
1421 // // position
1422 // node_pt(l)->x_gen(k,i) = old_var-fd_step;
1423
1424 // //Calculate the new residuals at backward
1425 // position get_residuals(newres_minus);
1426
1427 // //Do central finite differences
1428 // for(unsigned m=0;m<n_dof;m++)
1429 // {
1430 // //Stick the entry into the Jacobian matrix
1431 // jacobian(m,local_unknown) =
1432 // (newres[m] - newres_minus[m])/(2.0*fd_step);
1433 // }
1434 // }
1435
1436 // Reset the (generalised) Eulerian nodal position
1437 *value_pt = old_var;
1438
1439 // Perform any auxialiary node updates
1440 local_node_pt->perform_auxiliary_node_update_fct();
1441
1442 // Reset any other dependent variables
1444 }
1445 }
1446 }
1447 }
1448 // Otherwise it's a hanging node
1449 else
1450 {
1451 // Find the local hanging object
1452 HangInfo* hang_info_pt = local_node_pt->hanging_pt();
1453 // Loop over the master nodes
1454 const unsigned n_master = hang_info_pt->nmaster();
1455 for (unsigned m = 0; m < n_master; m++)
1456 {
1457 // Get the pointer to the master node
1458 Node* const master_node_pt = hang_info_pt->master_node_pt(m);
1459
1460 // Get the local equation numbers for the master node
1461 DenseMatrix<int> Position_local_eqn_at_node =
1462 Local_position_hang_eqn[master_node_pt];
1463
1464 // Loop over position dofs
1465 for (unsigned k = 0; k < n_position_type; k++)
1466 {
1467 // Loop over dimension
1468 for (unsigned i = 0; i < nodal_dim; i++)
1469 {
1470 local_unknown = Position_local_eqn_at_node(k, i);
1471 // If the variable is free
1472 if (local_unknown >= 0)
1473 {
1474 // Store a pointer to the (generalised) Eulerian nodal position
1475 double* const value_pt = &(master_node_pt->x_gen(k, i));
1476
1477 // Save the old value of the (generalised) Eulerian nodal
1478 // position
1479 const double old_var = *value_pt;
1480
1481 // Increment the (generalised) Eulerian nodal position
1482 *value_pt += fd_step;
1483
1484 // Perform any auxialiary node updates
1485 master_node_pt->perform_auxiliary_node_update_fct();
1486
1487 // Update any dependent variables
1489
1490 // Calculate the new residuals
1491 get_residuals(newres);
1492
1493 // if (use_first_order_fd)
1494 {
1495 // Do forward finite differences
1496 for (unsigned m = 0; m < n_dof; m++)
1497 {
1498 // Stick the entry into the Jacobian matrix
1499 jacobian(m, local_unknown) =
1500 (newres[m] - residuals[m]) / fd_step;
1501 }
1502 }
1503 // else
1504 // {
1505 // //Take backwards step for the (generalised)
1506 // Eulerian nodal
1507 // // position
1508 // master_node_pt->x_gen(k,i) = old_var-fd_step;
1509
1510 // //Calculate the new residuals at backward
1511 // position get_residuals(newres_minus);
1512
1513 // //Do central finite differences
1514 // for(unsigned m=0;m<n_dof;m++)
1515 // {
1516 // //Stick the entry into the Jacobian matrix
1517 // jacobian(m,local_unknown) =
1518 // (newres[m] -
1519 // newres_minus[m])/(2.0*fd_step);
1520 // }
1521 // }
1522
1523 // Reset the (generalised) Eulerian nodal position
1524 *value_pt = old_var;
1525
1526 // Perform any auxialiary node updates
1527 master_node_pt->perform_auxiliary_node_update_fct();
1528
1529 // Reset any other dependent variables
1531 }
1532 }
1533 }
1534 }
1535 } // End of hanging node case
1536
1537 } // End of loop over nodes
1538
1539 // End of finite difference loop
1540 // Final reset of any dependent data
1542 }
1543
1544} // namespace oomph
cstr elem_len * i
Definition: cfortran.h:603
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
Definition: shape.h:278
A class that represents a collection of data; each Data object may contain many different individual ...
Definition: nodes.h:86
static long Is_pinned
Static "Magic number" used in place of the equation number to indicate that the value is pinned.
Definition: nodes.h:183
unsigned nvalue() const
Return number of values stored in data object (incl pinned ones).
Definition: nodes.h:483
long & eqn_number(const unsigned &i)
Return the equation number of the i-th stored variable.
Definition: nodes.h:367
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:324
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:1709
unsigned nnodal_position_type() const
Return the number of coordinate types that the element requires to interpolate the geometry between t...
Definition: elements.h:2463
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:1718
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:1750
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:1432
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition: elements.h:2611
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:2349
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
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:1723
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:1714
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
unsigned nodal_dimension() const
Return the required Eulerian dimension of the nodes in this element.
Definition: elements.h:2484
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:1483
static double Default_fd_jacobian_step
Double used for the default finite difference step in elemental jacobian calculations.
Definition: elements.h:1198
unsigned ndof() const
Return the number of equations/dofs in the element.
Definition: elements.h:835
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:704
Data *& internal_data_pt(const unsigned &i)
Return a pointer to i-th internal data object.
Definition: elements.h:622
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:726
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:980
unsigned ninternal_data() const
Return the number of internal data objects.
Definition: elements.h:823
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:156
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:231
Class that contains data for hanging nodes.
Definition: nodes.h:742
Node *const & master_node_pt(const unsigned &i) const
Return a pointer to the i-th master node.
Definition: nodes.h:791
unsigned nmaster() const
Return the number of master nodes.
Definition: nodes.h:785
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
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:1615
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
Definition: nodes.h:1054
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:1060
bool is_hanging() const
Test whether the node is geometrically hanging.
Definition: nodes.h:1285
double & x_gen(const unsigned &k, const unsigned &i)
Reference to the generalised position x(k,i). ‘Type’: k; Coordinate direction: i.
Definition: nodes.h:1126
HangInfo *const & hanging_pt() const
Return pointer to hanging node data (this refers to the geometric hanging node status) (const version...
Definition: nodes.h:1228
An OomphLibError object which should be thrown when an run-time error is encountered....
////////////////////////////////////////////////////////////////// //////////////////////////////////...
Definition: matrices.h:1370
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...
static double Max_integrity_tolerance
Max. allowed discrepancy in element integrity check.
void assign_hanging_local_eqn_numbers(const bool &store_local_dof_pt)
Assign the local equation numbers for hanging node variables.
virtual void deactivate_element()
Final operations that must be performed when the element is no longer active in the mesh,...
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...
virtual unsigned ncont_interpolated_values() const =0
Number of continuously interpolated values. Note: We assume that they are located at the beginning of...
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.
void get_dresidual_dnodal_coordinates(RankThreeTensor< double > &dresidual_dnodal_coordinates)
Compute derivatives of elemental residual vector with respect to nodal coordinates....
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...
void assemble_local_to_eulerian_jacobian(const DShape &dpsids, DenseMatrix< double > &jacobian) const
Assemble the jacobian matrix for the mapping from local to Eulerian coordinates, given the derivative...
std::map< Node *, int > * Local_hang_eqn
Storage for local equation numbers of hanging node variables (values stored at master nodes)....
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 ...
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 identify_field_data_for_interactions(std::set< std::pair< Data *, unsigned > > &paired_field_data)
The purpose of this function is to identify all possible Data that can affect the fields interpolated...
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...
virtual ~RefineableElement()
Destructor, delete the allocated storage for the hanging equations.
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...
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...
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...
unsigned ngeom_data() const
The number of geometric data affecting a RefineableSolidFiniteElement is the positional Data of all n...
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...
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....
std::map< Node *, DenseMatrix< int > > Local_position_hang_eqn
Storage for local equation numbers of hanging node variables associated with nodal positions....
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...
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,...
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:3912
unsigned nnodal_lagrangian_type() const
Return the number of types of (generalised) nodal Lagrangian coordinates required to interpolate the ...
Definition: elements.h:3785
virtual void update_before_solid_position_fd()
Function that is called before the finite differencing of any solid position data....
Definition: elements.h:4240
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:4250
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:4137
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:4255
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:4245
A Class for nodes that deform elastically (i.e. position is an unknown in the problem)....
Definition: nodes.h:1686
A slight extension to the standard template vector class so that we can include "graceful" array rang...
Definition: Vector.h:58
//////////////////////////////////////////////////////////////////// ////////////////////////////////...