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