generalised_newtonian_refineable_navier_stokes_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// Header file for refineable 2D quad Navier Stokes elements
27
28#ifndef OOMPH_GENERALISED_NEWTONIAN_REFINEABLE_NAVIER_STOKES_ELEMENTS_HEADER
29#define OOMPH_GENERALISED_NEWTONIAN_REFINEABLE_NAVIER_STOKES_ELEMENTS_HEADER
30
31// Config header generated by autoconfig
32#ifdef HAVE_CONFIG_H
33#include <oomph-lib-config.h>
34#endif
35
36// Oomph-lib headers
37#include "../generic/refineable_quad_element.h"
38#include "../generic/refineable_brick_element.h"
39#include "../generic/hp_refineable_elements.h"
40#include "../generic/error_estimator.h"
42
43namespace oomph
44{
45 /// ////////////////////////////////////////////////////////////////////
46 /// ////////////////////////////////////////////////////////////////////
47 /// ////////////////////////////////////////////////////////////////////
48
49
50 //======================================================================
51 /// Refineable version of the Navier--Stokes equations
52 ///
53 ///
54 //======================================================================
55 template<unsigned DIM>
58 public virtual RefineableElement,
59 public virtual ElementWithZ2ErrorEstimator
60 {
61 protected:
62 /// Unpin all pressure dofs in the element
64
65 /// Pin unused nodal pressure dofs (empty by default, because
66 /// by default pressure dofs are not associated with nodes)
68
69 public:
70 /// Constructor
75 {
76 }
77
78
79 /// Loop over all elements in Vector (which typically contains
80 /// all the elements in a fluid mesh) and pin the nodal pressure degrees
81 /// of freedom that are not being used. Function uses
82 /// the member function
83 /// - \c RefineableNavierStokesEquations::
84 /// pin_elemental_redundant_nodal_pressure_dofs()
85 /// .
86 /// which is empty by default and should be implemented for
87 /// elements with nodal pressure degrees of freedom
88 /// (e.g. for refineable Taylor-Hood.)
90 const Vector<GeneralisedElement*>& element_pt)
91 {
92 // Loop over all elements and call the function that pins their
93 // unused nodal pressure data
94 unsigned n_element = element_pt.size();
95 for (unsigned e = 0; e < n_element; e++)
96 {
98 element_pt[e])
100 }
101 }
102
103 /// Unpin all pressure dofs in elements listed in vector.
105 const Vector<GeneralisedElement*>& element_pt)
106 {
107 // Loop over all elements
108 unsigned n_element = element_pt.size();
109 for (unsigned e = 0; e < n_element; e++)
110 {
112 element_pt[e])
114 }
115 }
116
117 /// Pointer to n_p-th pressure node (Default: NULL,
118 /// indicating that pressure is not based on nodal interpolation).
119 virtual Node* pressure_node_pt(const unsigned& n_p)
120 {
121 return NULL;
122 }
123
124 /// Compute the diagonal of the velocity/pressure mass matrices.
125 /// If which one=0, both are computed, otherwise only the pressure
126 /// (which_one=1) or the velocity mass matrix (which_one=2 -- the
127 /// LSC version of the preconditioner only needs that one)
129 Vector<double>& press_mass_diag,
130 Vector<double>& veloc_mass_diag,
131 const unsigned& which_one = 0);
132
133
134 /// Number of 'flux' terms for Z2 error estimation
136 {
137 // DIM diagonal strain rates, DIM(DIM -1) /2 off diagonal rates
138 return DIM + (DIM * (DIM - 1)) / 2;
139 }
140
141 /// Get 'flux' for Z2 error recovery: Upper triangular entries
142 /// in strain rate tensor.
144 {
145#ifdef PARANOID
146 unsigned num_entries = DIM + (DIM * (DIM - 1)) / 2;
147 if (flux.size() < num_entries)
148 {
149 std::ostringstream error_message;
150 error_message << "The flux vector has the wrong number of entries, "
151 << flux.size() << ", whereas it should be at least "
152 << num_entries << std::endl;
153 throw OomphLibError(error_message.str(),
154 OOMPH_CURRENT_FUNCTION,
155 OOMPH_EXCEPTION_LOCATION);
156 }
157#endif
158
159 // Get strain rate matrix
160 DenseMatrix<double> strainrate(DIM);
161 this->strain_rate(s, strainrate);
162
163 // Pack into flux Vector
164 unsigned icount = 0;
165
166 // Start with diagonal terms
167 for (unsigned i = 0; i < DIM; i++)
168 {
169 flux[icount] = strainrate(i, i);
170 icount++;
171 }
172
173 // Off diagonals row by row
174 for (unsigned i = 0; i < DIM; i++)
175 {
176 for (unsigned j = i + 1; j < DIM; j++)
177 {
178 flux[icount] = strainrate(i, j);
179 icount++;
180 }
181 }
182 }
183
184 /// Further build, pass the pointers down to the sons
186 {
187 // Find the father element
189 DIM>* cast_father_element_pt =
191 this->father_element_pt());
192
193 // Set the viscosity ratio pointer
194 this->Viscosity_Ratio_pt = cast_father_element_pt->viscosity_ratio_pt();
195 // Set the density ratio pointer
196 this->Density_Ratio_pt = cast_father_element_pt->density_ratio_pt();
197 // Set pointer to global Reynolds number
198 this->Re_pt = cast_father_element_pt->re_pt();
199 // Set pointer to global Reynolds number x Strouhal number (=Womersley)
200 this->ReSt_pt = cast_father_element_pt->re_st_pt();
201 // Set pointer to global Reynolds number x inverse Froude number
202 this->ReInvFr_pt = cast_father_element_pt->re_invfr_pt();
203 // Set pointer to global gravity Vector
204 this->G_pt = cast_father_element_pt->g_pt();
205
206 // Set pointer to body force function
207 this->Body_force_fct_pt = cast_father_element_pt->body_force_fct_pt();
208
209 // Set pointer to volumetric source function
210 this->Source_fct_pt = cast_father_element_pt->source_fct_pt();
211
212 // Set the ALE flag
213 this->ALE_is_disabled = cast_father_element_pt->ALE_is_disabled;
214 }
215
216
217 /// Compute the derivatives of the i-th component of
218 /// velocity at point s with respect
219 /// to all data that can affect its value. In addition, return the global
220 /// equation numbers corresponding to the data.
221 /// Overload the non-refineable version to take account of hanging node
222 /// information
224 const unsigned& i,
225 Vector<double>& du_ddata,
226 Vector<unsigned>& global_eqn_number)
227 {
228 // Find number of nodes
229 unsigned n_node = this->nnode();
230 // Local shape function
231 Shape psi(n_node);
232 // Find values of shape function at the given local coordinate
233 this->shape(s, psi);
234
235 // Find the index at which the velocity component is stored
236 const unsigned u_nodal_index = this->u_index_nst(i);
237
238 // Storage for hang info pointer
239 HangInfo* hang_info_pt = 0;
240 // Storage for global equation
241 int global_eqn = 0;
242
243 // Find the number of dofs associated with interpolated u
244 unsigned n_u_dof = 0;
245 for (unsigned l = 0; l < n_node; l++)
246 {
247 unsigned n_master = 1;
248
249 // Local bool (is the node hanging)
250 bool is_node_hanging = this->node_pt(l)->is_hanging();
251
252 // If the node is hanging, get the number of master nodes
253 if (is_node_hanging)
254 {
255 hang_info_pt = this->node_pt(l)->hanging_pt();
256 n_master = hang_info_pt->nmaster();
257 }
258 // Otherwise there is just one master node, the node itself
259 else
260 {
261 n_master = 1;
262 }
263
264 // Loop over the master nodes
265 for (unsigned m = 0; m < n_master; m++)
266 {
267 // Get the equation number
268 if (is_node_hanging)
269 {
270 // Get the equation number from the master node
271 global_eqn =
272 hang_info_pt->master_node_pt(m)->eqn_number(u_nodal_index);
273 }
274 else
275 {
276 // Global equation number
277 global_eqn = this->node_pt(l)->eqn_number(u_nodal_index);
278 }
279
280 // If it's positive add to the count
281 if (global_eqn >= 0)
282 {
283 ++n_u_dof;
284 }
285 }
286 }
287
288 // Now resize the storage schemes
289 du_ddata.resize(n_u_dof, 0.0);
290 global_eqn_number.resize(n_u_dof, 0);
291
292 // Loop over th nodes again and set the derivatives
293 unsigned count = 0;
294 // Loop over the local nodes and sum
295 for (unsigned l = 0; l < n_node; l++)
296 {
297 unsigned n_master = 1;
298 double hang_weight = 1.0;
299
300 // Local bool (is the node hanging)
301 bool is_node_hanging = this->node_pt(l)->is_hanging();
302
303 // If the node is hanging, get the number of master nodes
304 if (is_node_hanging)
305 {
306 hang_info_pt = this->node_pt(l)->hanging_pt();
307 n_master = hang_info_pt->nmaster();
308 }
309 // Otherwise there is just one master node, the node itself
310 else
311 {
312 n_master = 1;
313 }
314
315 // Loop over the master nodes
316 for (unsigned m = 0; m < n_master; m++)
317 {
318 // If the node is hanging get weight from master node
319 if (is_node_hanging)
320 {
321 // Get the hang weight from the master node
322 hang_weight = hang_info_pt->master_weight(m);
323 }
324 else
325 {
326 // Node contributes with full weight
327 hang_weight = 1.0;
328 }
329
330 // Get the equation number
331 if (is_node_hanging)
332 {
333 // Get the equation number from the master node
334 global_eqn =
335 hang_info_pt->master_node_pt(m)->eqn_number(u_nodal_index);
336 }
337 else
338 {
339 // Local equation number
340 global_eqn = this->node_pt(l)->eqn_number(u_nodal_index);
341 }
342
343 if (global_eqn >= 0)
344 {
345 // Set the global equation number
346 global_eqn_number[count] = global_eqn;
347 // Set the derivative with respect to the unknown
348 du_ddata[count] = psi[l] * hang_weight;
349 // Increase the counter
350 ++count;
351 }
352 }
353 }
354 }
355
356
357 protected:
358 /// Add element's contribution to elemental residual vector and/or
359 /// Jacobian matrix
360 /// flag=1: compute both
361 /// flag=0: compute only residual vector
363 Vector<double>& residuals,
364 DenseMatrix<double>& jacobian,
365 DenseMatrix<double>& mass_matrix,
366 unsigned flag);
367
368 /// Compute derivatives of elemental residual vector with respect
369 /// to nodal coordinates. Overwrites default implementation in
370 /// FiniteElement base class.
371 /// dresidual_dnodal_coordinates(l,i,j) = d res(l) / dX_{ij}
373 RankThreeTensor<double>& dresidual_dnodal_coordinates);
374 };
375
376
377 //======================================================================
378 /// Refineable version of Taylor Hood elements. These classes
379 /// can be written in total generality.
380 //======================================================================
381 template<unsigned DIM>
385 public virtual RefineableQElement<DIM>
386 {
387 private:
388 /// Unpin all pressure dofs
390 {
391 // find the index at which the pressure is stored
392 int p_index = this->p_nodal_index_nst();
393 unsigned n_node = this->nnode();
394 // loop over nodes
395 for (unsigned n = 0; n < n_node; n++)
396 {
397 this->node_pt(n)->unpin(p_index);
398 }
399 }
400
401 /// Pin all nodal pressure dofs that are not required
403 {
404 // Find the pressure index
405 int p_index = this->p_nodal_index_nst();
406 // Loop over all nodes
407 unsigned n_node = this->nnode();
408 // loop over all nodes and pin all the nodal pressures
409 for (unsigned n = 0; n < n_node; n++)
410 {
411 this->node_pt(n)->pin(p_index);
412 }
413
414 // Loop over all actual pressure nodes and unpin if they're not hanging
415 unsigned n_pres = this->npres_nst();
416 for (unsigned l = 0; l < n_pres; l++)
417 {
418 Node* nod_pt = this->pressure_node_pt(l);
419 if (!nod_pt->is_hanging(p_index))
420 {
421 nod_pt->unpin(p_index);
422 }
423 }
424 }
425
426 public:
427 /// Constructor
431 RefineableQElement<DIM>(),
433 {
434 }
435
436 /// Number of values required at local node n. In order to simplify
437 /// matters, we allocate storage for pressure variables at all the nodes
438 /// and then pin those that are not used.
439 unsigned required_nvalue(const unsigned& n) const
440 {
441 return DIM + 1;
442 }
443
444 /// Number of continuously interpolated values: (DIM velocities + 1
445 /// pressure)
447 {
448 return DIM + 1;
449 }
450
451 /// Rebuild from sons: empty
452 void rebuild_from_sons(Mesh*& mesh_pt) {}
453
454 /// Order of recovery shape functions for Z2 error estimation:
455 /// Same order as shape functions.
457 {
458 return 2;
459 }
460
461 /// Number of vertex nodes in the element
462 unsigned nvertex_node() const
463 {
465 }
466
467 /// Pointer to the j-th vertex node in the element
468 Node* vertex_node_pt(const unsigned& j) const
469 {
471 }
472
473 /// Get the function value u in Vector.
474 /// Note: Given the generality of the interface (this function
475 /// is usually called from black-box documentation or interpolation
476 /// routines), the values Vector sets its own size in here.
478 Vector<double>& values)
479 {
480 // Set size of Vector: u,v,p and initialise to zero
481 values.resize(DIM + 1, 0.0);
482
483 // Calculate velocities: values[0],...
484 for (unsigned i = 0; i < DIM; i++)
485 {
486 values[i] = this->interpolated_u_nst(s, i);
487 }
488
489 // Calculate pressure: values[DIM]
490 values[DIM] = this->interpolated_p_nst(s);
491 }
492
493 /// Get the function value u in Vector.
494 /// Note: Given the generality of the interface (this function
495 /// is usually called from black-box documentation or interpolation
496 /// routines), the values Vector sets its own size in here.
497 void get_interpolated_values(const unsigned& t,
498 const Vector<double>& s,
499 Vector<double>& values)
500 {
501 // Set size of Vector: u,v,p
502 values.resize(DIM + 1);
503
504 // Initialise
505 for (unsigned i = 0; i < DIM + 1; i++)
506 {
507 values[i] = 0.0;
508 }
509
510 // Find out how many nodes there are
511 unsigned n_node = this->nnode();
512
513 // Shape functions
514 Shape psif(n_node);
515 this->shape(s, psif);
516
517 // Calculate velocities: values[0],...
518 for (unsigned i = 0; i < DIM; i++)
519 {
520 // Get the index at which the i-th velocity is stored
521 unsigned u_nodal_index = this->u_index_nst(i);
522 for (unsigned l = 0; l < n_node; l++)
523 {
524 values[i] += this->nodal_value(t, l, u_nodal_index) * psif[l];
525 }
526 }
527
528 // Calculate pressure: values[DIM]
529 //(no history is carried in the pressure)
530 values[DIM] = this->interpolated_p_nst(s);
531 }
532
533
534 /// Perform additional hanging node procedures for variables
535 /// that are not interpolated by all nodes. The pressures are stored
536 /// at the p_nodal_index_nst-th location in each node
538 {
539 this->setup_hang_for_value(this->p_nodal_index_nst());
540 }
541
542 /// Pointer to n_p-th pressure node
543 Node* pressure_node_pt(const unsigned& n_p)
544 {
545 return this->node_pt(this->Pconv[n_p]);
546 }
547
548 /// The velocities are isoparametric and so the "nodes" interpolating
549 /// the velocities are the geometric nodes. The pressure "nodes" are a
550 /// subset of the nodes, so when value_id==DIM, the n-th pressure
551 /// node is returned.
552 Node* interpolating_node_pt(const unsigned& n, const int& value_id)
553
554 {
555 // The only different nodes are the pressure nodes
556 if (value_id == DIM)
557 {
558 return this->pressure_node_pt(n);
559 }
560 // The other variables are interpolated via the usual nodes
561 else
562 {
563 return this->node_pt(n);
564 }
565 }
566
567 /// The pressure nodes are the corner nodes, so when n_value==DIM,
568 /// the fraction is the same as the 1d node number, 0 or 1.
570 const unsigned& i,
571 const int& value_id)
572 {
573 if (value_id == DIM)
574 {
575 // The pressure nodes are just located on the boundaries at 0 or 1
576 return double(n1d);
577 }
578 // Otherwise the velocity nodes are the same as the geometric ones
579 else
580 {
581 return this->local_one_d_fraction_of_node(n1d, i);
582 }
583 }
584
585 /// The velocity nodes are the same as the geometric nodes. The
586 /// pressure nodes must be calculated by using the same methods as
587 /// the geometric nodes, but by recalling that there are only two pressure
588 /// nodes per edge.
590 const int& value_id)
591 {
592 // If we are calculating pressure nodes
593 if (value_id == DIM)
594 {
595 // Storage for the index of the pressure node
596 unsigned total_index = 0;
597 // The number of nodes along each 1d edge is 2.
598 unsigned NNODE_1D = 2;
599 // Storage for the index along each boundary
600 Vector<int> index(DIM);
601 // Loop over the coordinates
602 for (unsigned i = 0; i < DIM; i++)
603 {
604 // If we are at the lower limit, the index is zero
605 if (s[i] == -1.0)
606 {
607 index[i] = 0;
608 }
609 // If we are at the upper limit, the index is the number of nodes
610 // minus 1
611 else if (s[i] == 1.0)
612 {
613 index[i] = NNODE_1D - 1;
614 }
615 // Otherwise, we have to calculate the index in general
616 else
617 {
618 // For uniformly spaced nodes the 0th node number would be
619 double float_index = 0.5 * (1.0 + s[i]) * (NNODE_1D - 1);
620 index[i] = int(float_index);
621 // What is the excess. This should be safe because the
622 // taking the integer part rounds down
623 double excess = float_index - index[i];
624 // If the excess is bigger than our tolerance there is no node,
625 // return null
628 {
629 return 0;
630 }
631 }
632 /// Construct the general pressure index from the components.
633 total_index +=
634 index[i] * static_cast<unsigned>(pow(static_cast<float>(NNODE_1D),
635 static_cast<int>(i)));
636 }
637 // If we've got here we have a node, so let's return a pointer to it
638 return this->pressure_node_pt(total_index);
639 }
640 // Otherwise velocity nodes are the same as pressure nodes
641 else
642 {
643 return this->get_node_at_local_coordinate(s);
644 }
645 }
646
647
648 /// The number of 1d pressure nodes is 2, the number of 1d velocity
649 /// nodes is the same as the number of 1d geometric nodes.
650 unsigned ninterpolating_node_1d(const int& value_id)
651 {
652 if (value_id == DIM)
653 {
654 return 2;
655 }
656 else
657 {
658 return this->nnode_1d();
659 }
660 }
661
662 /// The number of pressure nodes is 2^DIM. The number of
663 /// velocity nodes is the same as the number of geometric nodes.
664 unsigned ninterpolating_node(const int& value_id)
665 {
666 if (value_id == DIM)
667 {
668 return static_cast<unsigned>(pow(2.0, static_cast<int>(DIM)));
669 }
670 else
671 {
672 return this->nnode();
673 }
674 }
675
676 /// The basis interpolating the pressure is given by pshape().
677 /// / The basis interpolating the velocity is shape().
679 Shape& psi,
680 const int& value_id) const
681 {
682 if (value_id == DIM)
683 {
684 return this->pshape_nst(s, psi);
685 }
686 else
687 {
688 return this->shape(s, psi);
689 }
690 }
691
692
693 /// Add to the set \c paired_load_data pairs containing
694 /// - the pointer to a Data object
695 /// and
696 /// - the index of the value in that Data object
697 /// .
698 /// for all values (pressures, velocities) that affect the
699 /// load computed in the \c get_load(...) function.
700 /// (Overloads non-refineable version and takes hanging nodes
701 /// into account)
703 std::set<std::pair<Data*, unsigned>>& paired_load_data)
704 {
705 // Get the nodal indices at which the velocities are stored
706 unsigned u_index[DIM];
707 for (unsigned i = 0; i < DIM; i++)
708 {
709 u_index[i] = this->u_index_nst(i);
710 }
711
712 // Loop over the nodes
713 unsigned n_node = this->nnode();
714 for (unsigned n = 0; n < n_node; n++)
715 {
716 // Pointer to current node
717 Node* nod_pt = this->node_pt(n);
718
719 // Check if it's hanging:
720 if (nod_pt->is_hanging())
721 {
722 // It's hanging -- get number of master nodes
723 unsigned nmaster = nod_pt->hanging_pt()->nmaster();
724
725 // Loop over masters
726 for (unsigned j = 0; j < nmaster; j++)
727 {
728 Node* master_nod_pt = nod_pt->hanging_pt()->master_node_pt(j);
729
730 // Loop over the velocity components and add pointer to their data
731 // and indices to the vectors
732 for (unsigned i = 0; i < DIM; i++)
733 {
734 paired_load_data.insert(
735 std::make_pair(master_nod_pt, u_index[i]));
736 }
737 }
738 }
739 // Not hanging
740 else
741 {
742 // Loop over the velocity components and add pointer to their data
743 // and indices to the vectors
744 for (unsigned i = 0; i < DIM; i++)
745 {
746 paired_load_data.insert(
747 std::make_pair(this->node_pt(n), u_index[i]));
748 }
749 }
750 }
751
752 // Get the nodal index at which the pressure is stored
753 int p_index = this->p_nodal_index_nst();
754
755 // Loop over the pressure data
756 unsigned n_pres = this->npres_nst();
757 for (unsigned l = 0; l < n_pres; l++)
758 {
759 // Get the pointer to the nodal pressure
760 Node* pres_node_pt = this->pressure_node_pt(l);
761 // Check if the pressure dof is hanging
762 if (pres_node_pt->is_hanging(p_index))
763 {
764 // Get the pointer to the hang info object
765 // (pressure is stored as p_index--th nodal dof).
766 HangInfo* hang_info_pt = pres_node_pt->hanging_pt(p_index);
767
768 // Get number of pressure master nodes (pressure is stored
769 unsigned nmaster = hang_info_pt->nmaster();
770
771 // Loop over pressure master nodes
772 for (unsigned m = 0; m < nmaster; m++)
773 {
774 // The p_index-th entry in each nodal data is the pressure, which
775 // affects the traction
776 paired_load_data.insert(
777 std::make_pair(hang_info_pt->master_node_pt(m), p_index));
778 }
779 }
780 // It's not hanging
781 else
782 {
783 // The p_index-th entry in each nodal data is the pressure, which
784 // affects the traction
785 paired_load_data.insert(std::make_pair(pres_node_pt, p_index));
786 }
787 }
788 }
789 };
790
791
792 //=======================================================================
793 /// Face geometry of the RefineableQTaylorHoodElements is the
794 /// same as the Face geometry of the QTaylorHoodElements.
795 //=======================================================================
796 template<unsigned DIM>
798 : public virtual FaceGeometry<GeneralisedNewtonianQTaylorHoodElement<DIM>>
799 {
800 public:
802 {
803 }
804 };
805
806
807 //=======================================================================
808 /// Face geometry of the face geometry of
809 /// the RefineableQTaylorHoodElements is the
810 /// same as the Face geometry of the Face geometry of QTaylorHoodElements.
811 //=======================================================================
812 template<unsigned DIM>
815 : public virtual FaceGeometry<
816 FaceGeometry<GeneralisedNewtonianQTaylorHoodElement<DIM>>>
817 {
818 public:
820 : FaceGeometry<
822 {
823 }
824 };
825
826
827 /// ////////////////////////////////////////////////////////////////////////
828 /// ////////////////////////////////////////////////////////////////////////
829 /// ////////////////////////////////////////////////////////////////////////
830
831
832 //======================================================================
833 /// Refineable version of Crouzeix Raviart elements. Generic class definitions
834 //======================================================================
835 template<unsigned DIM>
839 public virtual RefineableQElement<DIM>
840 {
841 private:
842 /// Unpin all internal pressure dofs
844 {
845 unsigned n_pres = this->npres_nst();
846 // loop over pressure dofs and unpin them
847 for (unsigned l = 0; l < n_pres; l++)
848 {
850 }
851 }
852
853 public:
854 /// Constructor
858 RefineableQElement<DIM>(),
860 {
861 }
862
863
864 /// Broken copy constructor
867 delete;
868
869 /// Broken assignment operator
870 // Commented out broken assignment operator because this can lead to a
871 // conflict warning when used in the virtual inheritence hierarchy.
872 // Essentially the compiler doesn't realise that two separate
873 // implementations of the broken function are the same and so, quite
874 // rightly, it shouts.
875 /*void operator=(const
876 RefineableGeneralisedNewtonianQCrouzeixRaviartElement<DIM>&)
877 = delete;*/
878
879 /// Number of continuously interpolated values: DIM (velocities)
881 {
882 return DIM;
883 }
884
885 /// Rebuild from sons: Reconstruct pressure from the (merged) sons
886 /// This must be specialised for each dimension.
887 inline void rebuild_from_sons(Mesh*& mesh_pt);
888
889 /// Order of recovery shape functions for Z2 error estimation:
890 /// Same order as shape functions.
892 {
893 return 2;
894 }
895
896 /// Number of vertex nodes in the element
897 unsigned nvertex_node() const
898 {
900 }
901
902 /// Pointer to the j-th vertex node in the element
903 Node* vertex_node_pt(const unsigned& j) const
904 {
906 j);
907 }
908
909 /// Get the function value u in Vector.
910 /// Note: Given the generality of the interface (this function
911 /// is usually called from black-box documentation or interpolation
912 /// routines), the values Vector sets its own size in here.
914 Vector<double>& values)
915 {
916 // Set size of Vector: u,v,p and initialise to zero
917 values.resize(DIM, 0.0);
918
919 // Calculate velocities: values[0],...
920 for (unsigned i = 0; i < DIM; i++)
921 {
922 values[i] = this->interpolated_u_nst(s, i);
923 }
924 }
925
926 /// Get all function values [u,v..,p] at previous timestep t
927 /// (t=0: present; t>0: previous timestep).
928 ///
929 /// Note: Given the generality of the interface (this function
930 /// is usually called from black-box documentation or interpolation
931 /// routines), the values Vector sets its own size in here.
932 ///
933 /// Note: No pressure history is kept, so pressure is always
934 /// the current value.
935 void get_interpolated_values(const unsigned& t,
936 const Vector<double>& s,
937 Vector<double>& values)
938 {
939 // Set size of Vector: u,v,p
940 values.resize(DIM);
941
942 // Initialise
943 for (unsigned i = 0; i < DIM; i++)
944 {
945 values[i] = 0.0;
946 }
947
948 // Find out how many nodes there are
949 unsigned n_node = this->nnode();
950
951 // Shape functions
952 Shape psif(n_node);
953 this->shape(s, psif);
954
955 // Calculate velocities: values[0],...
956 for (unsigned i = 0; i < DIM; i++)
957 {
958 // Get the nodal index at which the i-th velocity component is stored
959 unsigned u_nodal_index = this->u_index_nst(i);
960 for (unsigned l = 0; l < n_node; l++)
961 {
962 values[i] += this->nodal_value(t, l, u_nodal_index) * psif[l];
963 }
964 }
965 }
966
967 /// Perform additional hanging node procedures for variables
968 /// that are not interpolated by all nodes. Empty
970
971 /// Further build for Crouzeix_Raviart interpolates the internal
972 /// pressure dofs from father element: Make sure pressure values and
973 /// dp/ds agree between fathers and sons at the midpoints of the son
974 /// elements. This must be specialised for each dimension.
975 inline void further_build();
976
977
978 /// Add to the set \c paired_load_data pairs containing
979 /// - the pointer to a Data object
980 /// and
981 /// - the index of the value in that Data object
982 /// .
983 /// for all values (pressures, velocities) that affect the
984 /// load computed in the \c get_load(...) function.
985 /// (Overloads non-refineable version and takes hanging nodes
986 /// into account)
988 std::set<std::pair<Data*, unsigned>>& paired_load_data)
989 {
990 // Get the nodal indices at which the velocities are stored
991 unsigned u_index[DIM];
992 for (unsigned i = 0; i < DIM; i++)
993 {
994 u_index[i] = this->u_index_nst(i);
995 }
996
997 // Loop over the nodes
998 unsigned n_node = this->nnode();
999 for (unsigned n = 0; n < n_node; n++)
1000 {
1001 // Pointer to current node
1002 Node* nod_pt = this->node_pt(n);
1003
1004 // Check if it's hanging:
1005 if (nod_pt->is_hanging())
1006 {
1007 // It's hanging -- get number of master nodes
1008 unsigned nmaster = nod_pt->hanging_pt()->nmaster();
1009
1010 // Loop over masters
1011 for (unsigned j = 0; j < nmaster; j++)
1012 {
1013 Node* master_nod_pt = nod_pt->hanging_pt()->master_node_pt(j);
1014
1015 // Loop over the velocity components and add pointer to their data
1016 // and indices to the vectors
1017 for (unsigned i = 0; i < DIM; i++)
1018 {
1019 paired_load_data.insert(
1020 std::make_pair(master_nod_pt, u_index[i]));
1021 }
1022 }
1023 }
1024 // Not hanging
1025 else
1026 {
1027 // Loop over the velocity components and add pointer to their data
1028 // and indices to the vectors
1029 for (unsigned i = 0; i < DIM; i++)
1030 {
1031 paired_load_data.insert(
1032 std::make_pair(this->node_pt(n), u_index[i]));
1033 }
1034 }
1035 }
1036
1037
1038 // Loop over the pressure data (can't be hanging!)
1039 unsigned n_pres = this->npres_nst();
1040 for (unsigned l = 0; l < n_pres; l++)
1041 {
1042 // The entries in the internal data at P_nst_internal_index
1043 // are the pressures, which affect the traction
1044 paired_load_data.insert(std::make_pair(
1045 this->internal_data_pt(this->P_nst_internal_index), l));
1046 }
1047 }
1048 };
1049
1050
1051 //======================================================================
1052 /// p-refineable version of Crouzeix Raviart elements. Generic class
1053 /// definitions
1054 //======================================================================
1055 template<unsigned DIM>
1059 public virtual PRefineableQElement<DIM, 3>
1060 {
1061 private:
1062 /// Unpin all internal pressure dofs
1064 {
1065 unsigned n_pres = this->npres_nst();
1066 n_pres = this->internal_data_pt(this->P_nst_internal_index)->nvalue();
1067 // loop over pressure dofs and unpin them
1068 for (unsigned l = 0; l < n_pres; l++)
1069 {
1071 }
1072 }
1073
1074 public:
1075 /// Constructor
1079 PRefineableQElement<DIM, 3>(),
1081 {
1082 // Set the p-order
1083 this->p_order() = 3;
1084
1085 // Set integration scheme
1086 // (To avoid memory leaks in pre-build and p-refine where new
1087 // integration schemes are created)
1089
1090 // Resize pressure storage
1091 // (Constructor for QCrouzeixRaviartElement sets up DIM+1 pressure values)
1092 if (this->internal_data_pt(this->P_nst_internal_index)->nvalue() <=
1093 this->npres_nst())
1094 {
1096 ->resize(this->npres_nst());
1097 }
1098 else
1099 {
1100 Data* new_data_pt = new Data(this->npres_nst());
1101 delete this->internal_data_pt(this->P_nst_internal_index);
1102 this->internal_data_pt(this->P_nst_internal_index) = new_data_pt;
1103 }
1104 }
1105
1106 /// Destructor
1108 {
1109 delete this->integral_pt();
1110 }
1111
1112
1113 /// Broken copy constructor
1116 dummy) = delete;
1117
1118 /// Broken assignment operator
1119 /*void operator=(const
1120 PRefineableGeneralisedNewtonianQCrouzeixRaviartElement<DIM>&)
1121 = delete;*/
1122
1123 /// Return the i-th pressure value
1124 /// (Discontinous pressure interpolation -- no need to cater for hanging
1125 /// nodes).
1126 double p_nst(const unsigned& i) const
1127 {
1128 return this->internal_data_pt(this->P_nst_internal_index)->value(i);
1129 }
1130
1131 double p_nst(const unsigned& t, const unsigned& i) const
1132 {
1133 return this->internal_data_pt(this->P_nst_internal_index)->value(t, i);
1134 }
1135
1136 /// // Return number of pressure values
1137 unsigned npres_nst() const
1138 {
1139 return (this->p_order() - 2) * (this->p_order() - 2);
1140 }
1141
1142 /// Pin p_dof-th pressure dof and set it to value specified by p_value.
1143 void fix_pressure(const unsigned& p_dof, const double& p_value)
1144 {
1145 this->internal_data_pt(this->P_nst_internal_index)->pin(p_dof);
1147 ->set_value(p_dof, p_value);
1148 }
1149
1150 unsigned required_nvalue(const unsigned& n) const
1151 {
1152 return DIM;
1153 }
1154
1155 /// Number of continuously interpolated values: DIM (velocities)
1157 {
1158 return DIM;
1159 }
1160
1161 /// Rebuild from sons: Reconstruct pressure from the (merged) sons
1162 /// This must be specialised for each dimension.
1163 void rebuild_from_sons(Mesh*& mesh_pt)
1164 {
1165 // Do p-refineable version
1167 // Do Crouzeix-Raviart version
1168 // Need to reconstruct pressure manually!
1169 for (unsigned p = 0; p < npres_nst(); p++)
1170 {
1171 // BENFLAG: Set to zero for now -- don't do projection problem yet
1172 this->internal_data_pt(this->P_nst_internal_index)->set_value(p, 0.0);
1173 }
1174 }
1175
1176 /// Order of recovery shape functions for Z2 error estimation:
1177 /// - Same order as shape functions.
1178 // unsigned nrecovery_order()
1179 // {
1180 // if(this->nnode_1d() < 4) {return (this->nnode_1d()-1);}
1181 // else {return 3;}
1182 // }
1183 /// - Constant recovery order, since recovery order of the first element
1184 /// is used for the whole mesh.
1186 {
1187 return 3;
1188 }
1189
1190 /// Number of vertex nodes in the element
1191 unsigned nvertex_node() const
1192 {
1194 }
1195
1196 /// Pointer to the j-th vertex node in the element
1197 Node* vertex_node_pt(const unsigned& j) const
1198 {
1200 j);
1201 }
1202
1203 /// Velocity shape and test functions and their derivs
1204 /// w.r.t. to global coords at local coordinate s (taken from geometry)
1205 /// Return Jacobian of mapping between local and global coordinates.
1207 Shape& psi,
1208 DShape& dpsidx,
1209 Shape& test,
1210 DShape& dtestdx) const;
1211
1212 /// Velocity shape and test functions and their derivs
1213 /// w.r.t. to global coords at ipt-th integation point (taken from geometry)
1214 /// Return Jacobian of mapping between local and global coordinates.
1215 inline double dshape_and_dtest_eulerian_at_knot_nst(const unsigned& ipt,
1216 Shape& psi,
1217 DShape& dpsidx,
1218 Shape& test,
1219 DShape& dtestdx) const;
1220
1221 /// Pressure shape functions at local coordinate s
1222 inline void pshape_nst(const Vector<double>& s, Shape& psi) const;
1223
1224 /// Pressure shape and test functions at local coordinte s
1225 inline void pshape_nst(const Vector<double>& s,
1226 Shape& psi,
1227 Shape& test) const;
1228
1229 /// Get the function value u in Vector.
1230 /// Note: Given the generality of the interface (this function
1231 /// is usually called from black-box documentation or interpolation
1232 /// routines), the values Vector sets its own size in here.
1234 Vector<double>& values)
1235 {
1236 // Set size of Vector: u,v,p and initialise to zero
1237 values.resize(DIM, 0.0);
1238
1239 // Calculate velocities: values[0],...
1240 for (unsigned i = 0; i < DIM; i++)
1241 {
1242 values[i] = this->interpolated_u_nst(s, i);
1243 }
1244 }
1245
1246 /// Get all function values [u,v..,p] at previous timestep t
1247 /// (t=0: present; t>0: previous timestep).
1248 ///
1249 /// Note: Given the generality of the interface (this function
1250 /// is usually called from black-box documentation or interpolation
1251 /// routines), the values Vector sets its own size in here.
1252 ///
1253 /// Note: No pressure history is kept, so pressure is always
1254 /// the current value.
1255 void get_interpolated_values(const unsigned& t,
1256 const Vector<double>& s,
1257 Vector<double>& values)
1258 {
1259 // Set size of Vector: u,v,p
1260 values.resize(DIM);
1261
1262 // Initialise
1263 for (unsigned i = 0; i < DIM; i++)
1264 {
1265 values[i] = 0.0;
1266 }
1267
1268 // Find out how many nodes there are
1269 unsigned n_node = this->nnode();
1270
1271 // Shape functions
1272 Shape psif(n_node);
1273 this->shape(s, psif);
1274
1275 // Calculate velocities: values[0],...
1276 for (unsigned i = 0; i < DIM; i++)
1277 {
1278 // Get the nodal index at which the i-th velocity component is stored
1279 unsigned u_nodal_index = this->u_index_nst(i);
1280 for (unsigned l = 0; l < n_node; l++)
1281 {
1282 values[i] += this->nodal_value(t, l, u_nodal_index) * psif[l];
1283 }
1284 }
1285 }
1286
1287 /// Perform additional hanging node procedures for variables
1288 /// that are not interpolated by all nodes. Empty
1290
1291 /// Further build for Crouzeix_Raviart interpolates the internal
1292 /// pressure dofs from father element: Make sure pressure values and
1293 /// dp/ds agree between fathers and sons at the midpoints of the son
1294 /// elements. This must be specialised for each dimension.
1296 };
1297
1298
1299 //=======================================================================
1300 /// Face geometry of the RefineableQuadQCrouzeixRaviartElements
1301 //=======================================================================
1302 template<unsigned DIM>
1304 : public virtual FaceGeometry<
1305 GeneralisedNewtonianQCrouzeixRaviartElement<DIM>>
1306 {
1307 public:
1310 {
1311 }
1312 };
1313
1314 //======================================================================
1315 /// Face geometry of the face geometry of
1316 /// the RefineableQCrouzeixRaviartElements is the
1317 /// same as the Face geometry of the Face geometry of
1318 /// QCrouzeixRaviartElements.
1319 //=======================================================================
1320 template<unsigned DIM>
1323 : public virtual FaceGeometry<
1324 FaceGeometry<GeneralisedNewtonianQCrouzeixRaviartElement<DIM>>>
1325 {
1326 public:
1328 : FaceGeometry<
1330 {
1331 }
1332 };
1333
1334
1335 // Inline functions
1336
1337 //=====================================================================
1338 /// 2D Rebuild from sons: Reconstruct pressure from the (merged) sons
1339 //=====================================================================
1340 template<>
1341 inline void RefineableGeneralisedNewtonianQCrouzeixRaviartElement<
1342 2>::rebuild_from_sons(Mesh*& mesh_pt)
1343 {
1344 using namespace QuadTreeNames;
1345
1346 // Central pressure value:
1347 //-----------------------
1348
1349 // Use average of the sons central pressure values
1350 // Other options: Take average of the four (discontinuous)
1351 // pressure values at the father's midpoint]
1352
1353 double av_press = 0.0;
1354
1355 // Loop over the sons
1356 for (unsigned ison = 0; ison < 4; ison++)
1357 {
1358 // Add the sons midnode pressure
1359 // Note that we can assume that the pressure is stored at the same
1360 // location because these are EXACTLY the same type of elements
1361 av_press += quadtree_pt()
1362 ->son_pt(ison)
1363 ->object_pt()
1364 ->internal_data_pt(this->P_nst_internal_index)
1365 ->value(0);
1366 }
1367
1368 // Use the average
1369 internal_data_pt(this->P_nst_internal_index)->set_value(0, 0.25 * av_press);
1370
1371
1372 // Slope in s_0 direction
1373 //----------------------
1374
1375 // Use average of the 2 FD approximations based on the
1376 // elements central pressure values
1377 // [Other options: Take average of the four
1378 // pressure derivatives]
1379
1380 double slope1 = quadtree_pt()
1381 ->son_pt(SE)
1382 ->object_pt()
1383 ->internal_data_pt(this->P_nst_internal_index)
1384 ->value(0) -
1385 quadtree_pt()
1386 ->son_pt(SW)
1387 ->object_pt()
1388 ->internal_data_pt(this->P_nst_internal_index)
1389 ->value(0);
1390
1391 double slope2 = quadtree_pt()
1392 ->son_pt(NE)
1393 ->object_pt()
1394 ->internal_data_pt(this->P_nst_internal_index)
1395 ->value(0) -
1396 quadtree_pt()
1397 ->son_pt(NW)
1398 ->object_pt()
1399 ->internal_data_pt(this->P_nst_internal_index)
1400 ->value(0);
1401
1402
1403 // Use the average
1404 internal_data_pt(this->P_nst_internal_index)
1405 ->set_value(1, 0.5 * (slope1 + slope2));
1406
1407
1408 // Slope in s_1 direction
1409 //----------------------
1410
1411 // Use average of the 2 FD approximations based on the
1412 // elements central pressure values
1413 // [Other options: Take average of the four
1414 // pressure derivatives]
1415
1416 slope1 = quadtree_pt()
1417 ->son_pt(NE)
1418 ->object_pt()
1419 ->internal_data_pt(this->P_nst_internal_index)
1420 ->value(0) -
1421 quadtree_pt()
1422 ->son_pt(SE)
1423 ->object_pt()
1424 ->internal_data_pt(this->P_nst_internal_index)
1425 ->value(0);
1426
1427 slope2 = quadtree_pt()
1428 ->son_pt(NW)
1429 ->object_pt()
1430 ->internal_data_pt(this->P_nst_internal_index)
1431 ->value(0) -
1432 quadtree_pt()
1433 ->son_pt(SW)
1434 ->object_pt()
1435 ->internal_data_pt(this->P_nst_internal_index)
1436 ->value(0);
1437
1438
1439 // Use the average
1440 internal_data_pt(this->P_nst_internal_index)
1441 ->set_value(2, 0.5 * (slope1 + slope2));
1442 }
1443
1444
1445 //=================================================================
1446 /// 3D Rebuild from sons: Reconstruct pressure from the (merged) sons
1447 //=================================================================
1448 template<>
1450 3>::rebuild_from_sons(Mesh*& mesh_pt)
1451 {
1452 using namespace OcTreeNames;
1453
1454 // Central pressure value:
1455 //-----------------------
1456
1457 // Use average of the sons central pressure values
1458 // Other options: Take average of the four (discontinuous)
1459 // pressure values at the father's midpoint]
1460
1461 double av_press = 0.0;
1462
1463 // Loop over the sons
1464 for (unsigned ison = 0; ison < 8; ison++)
1465 {
1466 // Add the sons midnode pressure
1467 av_press += octree_pt()
1468 ->son_pt(ison)
1469 ->object_pt()
1470 ->internal_data_pt(this->P_nst_internal_index)
1471 ->value(0);
1472 }
1473
1474 // Use the average
1475 internal_data_pt(this->P_nst_internal_index)
1476 ->set_value(0, 0.125 * av_press);
1477
1478
1479 // Slope in s_0 direction
1480 //----------------------
1481
1482 // Use average of the 4 FD approximations based on the
1483 // elements central pressure values
1484 // [Other options: Take average of the four
1485 // pressure derivatives]
1486
1487 double slope1 = octree_pt()
1488 ->son_pt(RDF)
1489 ->object_pt()
1490 ->internal_data_pt(this->P_nst_internal_index)
1491 ->value(0) -
1492 octree_pt()
1493 ->son_pt(LDF)
1494 ->object_pt()
1495 ->internal_data_pt(this->P_nst_internal_index)
1496 ->value(0);
1497
1498 double slope2 = octree_pt()
1499 ->son_pt(RUF)
1500 ->object_pt()
1501 ->internal_data_pt(this->P_nst_internal_index)
1502 ->value(0) -
1503 octree_pt()
1504 ->son_pt(LUF)
1505 ->object_pt()
1506 ->internal_data_pt(this->P_nst_internal_index)
1507 ->value(0);
1508
1509 double slope3 = octree_pt()
1510 ->son_pt(RDB)
1511 ->object_pt()
1512 ->internal_data_pt(this->P_nst_internal_index)
1513 ->value(0) -
1514 octree_pt()
1515 ->son_pt(LDB)
1516 ->object_pt()
1517 ->internal_data_pt(this->P_nst_internal_index)
1518 ->value(0);
1519
1520 double slope4 = octree_pt()
1521 ->son_pt(RUB)
1522 ->object_pt()
1523 ->internal_data_pt(this->P_nst_internal_index)
1524 ->value(0) -
1525 octree_pt()
1526 ->son_pt(LUB)
1527 ->object_pt()
1528 ->internal_data_pt(this->P_nst_internal_index)
1529 ->value(0);
1530
1531
1532 // Use the average
1533 internal_data_pt(this->P_nst_internal_index)
1534 ->set_value(1, 0.25 * (slope1 + slope2 + slope3 + slope4));
1535
1536
1537 // Slope in s_1 direction
1538 //----------------------
1539
1540 // Use average of the 4 FD approximations based on the
1541 // elements central pressure values
1542 // [Other options: Take average of the four
1543 // pressure derivatives]
1544
1545 slope1 = octree_pt()
1546 ->son_pt(LUB)
1547 ->object_pt()
1548 ->internal_data_pt(this->P_nst_internal_index)
1549 ->value(0) -
1550 octree_pt()
1551 ->son_pt(LDB)
1552 ->object_pt()
1553 ->internal_data_pt(this->P_nst_internal_index)
1554 ->value(0);
1555
1556 slope2 = octree_pt()
1557 ->son_pt(RUB)
1558 ->object_pt()
1559 ->internal_data_pt(this->P_nst_internal_index)
1560 ->value(0) -
1561 octree_pt()
1562 ->son_pt(RDB)
1563 ->object_pt()
1564 ->internal_data_pt(this->P_nst_internal_index)
1565 ->value(0);
1566
1567 slope3 = octree_pt()
1568 ->son_pt(LUF)
1569 ->object_pt()
1570 ->internal_data_pt(this->P_nst_internal_index)
1571 ->value(0) -
1572 octree_pt()
1573 ->son_pt(LDF)
1574 ->object_pt()
1575 ->internal_data_pt(this->P_nst_internal_index)
1576 ->value(0);
1577
1578 slope4 = octree_pt()
1579 ->son_pt(RUF)
1580 ->object_pt()
1581 ->internal_data_pt(this->P_nst_internal_index)
1582 ->value(0) -
1583 octree_pt()
1584 ->son_pt(RDF)
1585 ->object_pt()
1586 ->internal_data_pt(this->P_nst_internal_index)
1587 ->value(0);
1588
1589
1590 // Use the average
1591 internal_data_pt(this->P_nst_internal_index)
1592 ->set_value(2, 0.25 * (slope1 + slope2 + slope3 + slope4));
1593
1594
1595 // Slope in s_2 direction
1596 //----------------------
1597
1598 // Use average of the 4 FD approximations based on the
1599 // elements central pressure values
1600 // [Other options: Take average of the four
1601 // pressure derivatives]
1602
1603 slope1 = octree_pt()
1604 ->son_pt(LUF)
1605 ->object_pt()
1606 ->internal_data_pt(this->P_nst_internal_index)
1607 ->value(0) -
1608 octree_pt()
1609 ->son_pt(LUB)
1610 ->object_pt()
1611 ->internal_data_pt(this->P_nst_internal_index)
1612 ->value(0);
1613
1614 slope2 = octree_pt()
1615 ->son_pt(RUF)
1616 ->object_pt()
1617 ->internal_data_pt(this->P_nst_internal_index)
1618 ->value(0) -
1619 octree_pt()
1620 ->son_pt(RUB)
1621 ->object_pt()
1622 ->internal_data_pt(this->P_nst_internal_index)
1623 ->value(0);
1624
1625 slope3 = octree_pt()
1626 ->son_pt(LDF)
1627 ->object_pt()
1628 ->internal_data_pt(this->P_nst_internal_index)
1629 ->value(0) -
1630 octree_pt()
1631 ->son_pt(LDB)
1632 ->object_pt()
1633 ->internal_data_pt(this->P_nst_internal_index)
1634 ->value(0);
1635
1636 slope4 = octree_pt()
1637 ->son_pt(RDF)
1638 ->object_pt()
1639 ->internal_data_pt(this->P_nst_internal_index)
1640 ->value(0) -
1641 octree_pt()
1642 ->son_pt(RDB)
1643 ->object_pt()
1644 ->internal_data_pt(this->P_nst_internal_index)
1645 ->value(0);
1646
1647 // Use the average
1648 internal_data_pt(this->P_nst_internal_index)
1649 ->set_value(3, 0.25 * (slope1 + slope2 + slope3 + slope4));
1650 }
1651
1652
1653 //======================================================================
1654 /// 2D Further build for Crouzeix_Raviart interpolates the internal
1655 /// pressure dofs from father element: Make sure pressure values and
1656 /// dp/ds agree between fathers and sons at the midpoints of the son
1657 /// elements.
1658 //======================================================================
1659 template<>
1661 2>::further_build()
1662 {
1663 // Call the generic further build
1665
1666 using namespace QuadTreeNames;
1667
1668 // What type of son am I? Ask my quadtree representation...
1669 int son_type = quadtree_pt()->son_type();
1670
1671 // Pointer to my father (in element impersonation)
1672 RefineableElement* father_el_pt = quadtree_pt()->father_pt()->object_pt();
1673
1674 Vector<double> s_father(2);
1675
1676 // Son midpoint is located at the following coordinates in father element:
1677
1678 // South west son
1679 if (son_type == SW)
1680 {
1681 s_father[0] = -0.5;
1682 s_father[1] = -0.5;
1683 }
1684 // South east son
1685 else if (son_type == SE)
1686 {
1687 s_father[0] = 0.5;
1688 s_father[1] = -0.5;
1689 }
1690 // North east son
1691 else if (son_type == NE)
1692 {
1693 s_father[0] = 0.5;
1694 s_father[1] = 0.5;
1695 }
1696
1697 // North west son
1698 else if (son_type == NW)
1699 {
1700 s_father[0] = -0.5;
1701 s_father[1] = 0.5;
1702 }
1703
1704 // Pressure value in father element
1706 cast_father_element_pt =
1708 father_el_pt);
1709
1710 double press = cast_father_element_pt->interpolated_p_nst(s_father);
1711
1712 // Pressure value gets copied straight into internal dof:
1713 internal_data_pt(this->P_nst_internal_index)->set_value(0, press);
1714
1715 // The slopes get copied from father
1716 for (unsigned i = 1; i < 3; i++)
1717 {
1718 double half_father_slope =
1719 0.5 *
1720 cast_father_element_pt->internal_data_pt(this->P_nst_internal_index)
1721 ->value(i);
1722 // Set the value in the son
1723 internal_data_pt(this->P_nst_internal_index)
1724 ->set_value(i, half_father_slope);
1725 }
1726 }
1727
1728
1729 //=======================================================================
1730 /// 3D Further build for Crouzeix_Raviart interpolates the internal
1731 /// pressure dofs from father element: Make sure pressure values and
1732 /// dp/ds agree between fathers and sons at the midpoints of the son
1733 /// elements.
1734 //=======================================================================
1735 template<>
1737 3>::further_build()
1738 {
1740
1741 using namespace OcTreeNames;
1742
1743 // What type of son am I? Ask my octree representation...
1744 int son_type = octree_pt()->son_type();
1745
1746 // Pointer to my father (in element impersonation)
1747 RefineableQElement<3>* father_el_pt = dynamic_cast<RefineableQElement<3>*>(
1748 octree_pt()->father_pt()->object_pt());
1749
1750 Vector<double> s_father(3);
1751
1752 // Son midpoint is located at the following coordinates in father element:
1753 for (unsigned i = 0; i < 3; i++)
1754 {
1755 s_father[i] = 0.5 * OcTree::Direction_to_vector[son_type][i];
1756 }
1757
1758 // Pressure value in father element
1760 cast_father_element_pt =
1762 father_el_pt);
1763
1764 double press = cast_father_element_pt->interpolated_p_nst(s_father);
1765
1766 // Pressure value gets copied straight into internal dof:
1767 internal_data_pt(this->P_nst_internal_index)->set_value(0, press);
1768
1769 // The slopes get copied from father
1770 for (unsigned i = 1; i < 4; i++)
1771 {
1772 double half_father_slope =
1773 0.5 *
1774 cast_father_element_pt->internal_data_pt(this->P_nst_internal_index)
1775 ->value(i);
1776 // Set the value
1777 internal_data_pt(this->P_nst_internal_index)
1778 ->set_value(i, half_father_slope);
1779 }
1780 }
1781
1782 //=======================================================================
1783 /// 2D
1784 /// Derivatives of the shape functions and test functions w.r.t. to global
1785 /// (Eulerian) coordinates. Return Jacobian of mapping between
1786 /// local and global coordinates.
1787 //=======================================================================
1788 template<>
1790 2>::dshape_and_dtest_eulerian_nst(const Vector<double>& s,
1791 Shape& psi,
1792 DShape& dpsidx,
1793 Shape& test,
1794 DShape& dtestdx) const
1795 {
1796 // Call the geometrical shape functions and derivatives
1797 double J = this->dshape_eulerian(s, psi, dpsidx);
1798
1799 // Loop over the test functions and derivatives and set them equal to the
1800 // shape functions
1801 for (unsigned i = 0; i < nnode_1d() * nnode_1d(); i++)
1802 {
1803 test[i] = psi[i];
1804 dtestdx(i, 0) = dpsidx(i, 0);
1805 dtestdx(i, 1) = dpsidx(i, 1);
1806 }
1807
1808 // Return the jacobian
1809 return J;
1810 }
1811
1812 //=======================================================================
1813 /// 2D
1814 /// Derivatives of the shape functions and test functions w.r.t. to global
1815 /// (Eulerian) coordinates. Return Jacobian of mapping between
1816 /// local and global coordinates.
1817 //=======================================================================
1818 template<>
1820 2>::dshape_and_dtest_eulerian_at_knot_nst(const unsigned& ipt,
1821 Shape& psi,
1822 DShape& dpsidx,
1823 Shape& test,
1824 DShape& dtestdx) const
1825 {
1826 // Call the geometrical shape functions and derivatives
1827 double J = this->dshape_eulerian_at_knot(ipt, psi, dpsidx);
1828
1829 // Loop over the test functions and derivatives and set them equal to the
1830 // shape functions
1831 for (unsigned i = 0; i < nnode_1d() * nnode_1d(); i++)
1832 {
1833 test[i] = psi[i];
1834 dtestdx(i, 0) = dpsidx(i, 0);
1835 dtestdx(i, 1) = dpsidx(i, 1);
1836 }
1837
1838 // Return the jacobian
1839 return J;
1840 }
1841
1842 //=======================================================================
1843 /// 3D
1844 /// Derivatives of the shape functions and test functions w.r.t. to global
1845 /// (Eulerian) coordinates. Return Jacobian of mapping between
1846 /// local and global coordinates.
1847 //=======================================================================
1848 template<>
1850 3>::dshape_and_dtest_eulerian_nst(const Vector<double>& s,
1851 Shape& psi,
1852 DShape& dpsidx,
1853 Shape& test,
1854 DShape& dtestdx) const
1855 {
1856 // Call the geometrical shape functions and derivatives
1857 double J = this->dshape_eulerian(s, psi, dpsidx);
1858
1859 // Loop over the test functions and derivatives and set them equal to the
1860 // shape functions
1861 for (unsigned i = 0; i < nnode_1d() * nnode_1d() * nnode_1d(); i++)
1862 {
1863 test[i] = psi[i];
1864 dtestdx(i, 0) = dpsidx(i, 0);
1865 dtestdx(i, 1) = dpsidx(i, 1);
1866 dtestdx(i, 2) = dpsidx(i, 2);
1867 }
1868
1869 // Return the jacobian
1870 return J;
1871 }
1872
1873 //=======================================================================
1874 /// 3D
1875 /// Derivatives of the shape functions and test functions w.r.t. to global
1876 /// (Eulerian) coordinates. Return Jacobian of mapping between
1877 /// local and global coordinates.
1878 //=======================================================================
1879 template<>
1881 3>::dshape_and_dtest_eulerian_at_knot_nst(const unsigned& ipt,
1882 Shape& psi,
1883 DShape& dpsidx,
1884 Shape& test,
1885 DShape& dtestdx) const
1886 {
1887 // Call the geometrical shape functions and derivatives
1888 double J = this->dshape_eulerian_at_knot(ipt, psi, dpsidx);
1889
1890 // Loop over the test functions and derivatives and set them equal to the
1891 // shape functions
1892 for (unsigned i = 0; i < nnode_1d() * nnode_1d() * nnode_1d(); i++)
1893 {
1894 test[i] = psi[i];
1895 dtestdx(i, 0) = dpsidx(i, 0);
1896 dtestdx(i, 1) = dpsidx(i, 1);
1897 dtestdx(i, 2) = dpsidx(i, 2);
1898 }
1899
1900 // Return the jacobian
1901 return J;
1902 }
1903
1904 //=======================================================================
1905 /// 2D :
1906 /// Pressure shape functions
1907 //=======================================================================
1908 template<>
1910 2>::pshape_nst(const Vector<double>& s, Shape& psi) const
1911 {
1912 unsigned npres = this->npres_nst();
1913 if (npres == 1)
1914 {
1915 psi[0] = 1.0;
1916 }
1917 else
1918 {
1919 // Get number of pressure modes
1920 unsigned npres_1d = (int)std::sqrt((double)npres);
1921
1922 // Local storage
1923 // Call the one-dimensional modal shape functions
1924 OneDimensionalModalShape psi1(npres_1d, s[0]);
1925 OneDimensionalModalShape psi2(npres_1d, s[1]);
1926
1927 // Now let's loop over the nodal points in the element
1928 // s1 is the "x" coordinate, s2 the "y"
1929 for (unsigned i = 0; i < npres_1d; i++)
1930 {
1931 for (unsigned j = 0; j < npres_1d; j++)
1932 {
1933 // Multiply the two 1D functions together to get the 2D function
1934 psi[i * npres_1d + j] = psi2[i] * psi1[j];
1935 }
1936 }
1937 }
1938 }
1939
1940 /// Define the pressure shape and test functions
1941 template<>
1943 2>::pshape_nst(const Vector<double>& s, Shape& psi, Shape& test) const
1944 {
1945 // Call the pressure shape functions
1946 pshape_nst(s, psi);
1947
1948 // Loop over the test functions and set them equal to the shape functions
1949 if (this->npres_nst() == 1)
1950 {
1951 test[0] = psi[0];
1952 }
1953 else
1954 {
1955 for (unsigned i = 0; i < this->npres_nst(); i++) test[i] = psi[i];
1956 }
1957 }
1958
1959 //=======================================================================
1960 /// 3D :
1961 /// Pressure shape functions
1962 //=======================================================================
1963 template<>
1965 3>::pshape_nst(const Vector<double>& s, Shape& psi) const
1966 {
1967 unsigned npres = this->npres_nst();
1968 if (npres == 1)
1969 {
1970 psi[0] = 1.0;
1971 }
1972 else
1973 {
1974 // Get number of pressure modes
1975 unsigned npres_1d = (int)std::sqrt((double)npres);
1976
1977 // Local storage
1978 // Call the one-dimensional modal shape functions
1979 OneDimensionalModalShape psi1(npres_1d, s[0]);
1980 OneDimensionalModalShape psi2(npres_1d, s[1]);
1981 OneDimensionalModalShape psi3(npres_1d, s[2]);
1982
1983 // Now let's loop over the nodal points in the element
1984 // s1 is the "x" coordinate, s2 the "y"
1985 for (unsigned i = 0; i < npres_1d; i++)
1986 {
1987 for (unsigned j = 0; j < npres_1d; j++)
1988 {
1989 for (unsigned k = 0; k < npres_1d; k++)
1990 {
1991 // Multiply the two 1D functions together to get the 2D function
1992 psi[i * npres_1d * npres_1d + j * npres_1d + k] =
1993 psi3[i] * psi2[j] * psi1[k];
1994 }
1995 }
1996 }
1997 }
1998 }
1999
2000 /// Define the pressure shape and test functions
2001 template<>
2003 3>::pshape_nst(const Vector<double>& s, Shape& psi, Shape& test) const
2004 {
2005 // Call the pressure shape functions
2006 pshape_nst(s, psi);
2007
2008 // Loop over the test functions and set them equal to the shape functions
2009 if (this->npres_nst() == 1)
2010 {
2011 test[0] = psi[0];
2012 }
2013 else
2014 {
2015 for (unsigned i = 0; i < this->npres_nst(); i++) test[i] = psi[i];
2016 }
2017 }
2018
2019} // namespace oomph
2020
2021#endif
e
Definition: cfortran.h:571
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
char t
Definition: cfortran.h:568
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
void pin(const unsigned &i)
Pin the i-th stored variable.
Definition: nodes.h:385
void unpin(const unsigned &i)
Unpin the i-th stored variable.
Definition: nodes.h:391
void set_value(const unsigned &i, const double &value_)
Set the i-th stored data value to specified value. The only reason that we require an explicit set fu...
Definition: nodes.h:271
unsigned nvalue() const
Return number of values stored in data object (incl pinned ones).
Definition: nodes.h:483
double value(const unsigned &i) const
Return i-th stored value. This function is not virtual so that it can be inlined. This means that if ...
Definition: nodes.h:293
long & eqn_number(const unsigned &i)
Return the equation number of the i-th stored variable.
Definition: nodes.h:367
virtual void resize(const unsigned &n_value)
Change (increase) the number of values that may be stored.
Definition: nodes.cc:1002
Base class for finite elements that can compute the quantities that are required for the Z2 error est...
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: elements.h:4998
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1963
double nodal_value(const unsigned &n, const unsigned &i) const
Return the i-th value stored at local node n. Produces suitably interpolated values for hanging nodes...
Definition: elements.h:2593
virtual Node * get_node_at_local_coordinate(const Vector< double > &s) const
If there is a node at this local coordinate, return the pointer to the node.
Definition: elements.cc:3882
virtual unsigned nvertex_node() const
Return the number of vertex nodes in this element. Broken virtual function in "pure" finite elements.
Definition: elements.h:2491
static const double Node_location_tolerance
Default value for the tolerance to be used when locating nodes via local coordinates.
Definition: elements.h:1374
virtual void shape(const Vector< double > &s, Shape &psi) const =0
Calculate the geometric shape functions at local coordinate s. This function must be overloaded for e...
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
virtual unsigned nnode_1d() const
Return the number of nodes along one edge of the element Default is to return zero — must be overload...
Definition: elements.h:2218
virtual void set_integration_scheme(Integral *const &integral_pt)
Set the spatial integration scheme.
Definition: elements.cc:3210
virtual Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element. Broken virtual function in "pure" finite elements.
Definition: elements.h:2500
virtual double local_one_d_fraction_of_node(const unsigned &n1d, const unsigned &i)
Get the local fraction of any node in the n-th position in a one dimensional expansion along the i-th...
Definition: elements.h:1858
Class for multidimensional Gauss Lobatto Legendre integration rules empty - just establishes template...
Definition: integral.h:1281
Data *& internal_data_pt(const unsigned &i)
Return a pointer to i-th internal data object.
Definition: elements.h:622
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
double * ReInvFr_pt
Pointer to global Reynolds number x inverse Froude number (= Bond number / Capillary number)
double * Viscosity_Ratio_pt
Pointer to the viscosity ratio (relative to the viscosity used in the definition of the Reynolds numb...
bool ALE_is_disabled
Boolean flag to indicate if ALE formulation is disabled when time-derivatives are computed....
void strain_rate(const Vector< double > &s, DenseMatrix< double > &strain_rate) const
Strain-rate tensor: 1/2 (du_i/dx_j + du_j/dx_i)
double * Density_Ratio_pt
Pointer to the density ratio (relative to the density used in the definition of the Reynolds number)
virtual unsigned u_index_nst(const unsigned &i) const
Return the index at which the i-th unknown velocity component is stored. The default value,...
void interpolated_u_nst(const Vector< double > &s, Vector< double > &veloc) const
Compute vector of FE interpolated velocity u at local coordinate s.
NavierStokesBodyForceFctPt Body_force_fct_pt
Pointer to body force function.
NavierStokesSourceFctPt Source_fct_pt
Pointer to volumetric source function.
double * ReSt_pt
Pointer to global Reynolds number x Strouhal number (=Womersley)
virtual double interpolated_p_nst(const Vector< double > &s) const
Return FE interpolated pressure at local coordinate s.
/////////////////////////////////////////////////////////////////////////// /////////////////////////...
unsigned P_nst_internal_index
Internal index that indicates at which internal data the pressure is stored.
/////////////////////////////////////////////////////////////////////////
virtual int p_nodal_index_nst() const
Set the value at which the pressure is stored in the nodes.
void pshape_nst(const Vector< double > &s, Shape &psi) const
Pressure shape functions at local coordinate s.
static const unsigned Pconv[]
Static array of ints to hold conversion from pressure node numbers to actual node numbers.
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
double const & master_weight(const unsigned &i) const
Return weight for dofs on i-th master node.
Definition: nodes.h:808
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
bool is_hanging() const
Test whether the node is geometrically hanging.
Definition: nodes.h:1285
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
static Vector< Vector< int > > Direction_to_vector
For each direction, i.e. a son_type (vertex), a face or an edge, this defines a vector that indicates...
Definition: octree.h:353
Non-templated class that returns modal hierachical shape functions based on Legendre polynomials.
Definition: shape.h:1349
An OomphLibError object which should be thrown when an run-time error is encountered....
p-refineable version of Crouzeix Raviart elements. Generic class definitions
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation:
void pshape_nst(const Vector< double > &s, Shape &psi) const
Pressure shape functions at local coordinate s.
void further_setup_hanging_nodes()
Perform additional hanging node procedures for variables that are not interpolated by all nodes....
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
double dshape_and_dtest_eulerian_at_knot_nst(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const
Velocity shape and test functions and their derivs w.r.t. to global coords at ipt-th integation point...
double p_nst(const unsigned &t, const unsigned &i) const
Return the i-th pressure value (Discontinous pressure interpolation – no need to cater for hanging no...
void rebuild_from_sons(Mesh *&mesh_pt)
Rebuild from sons: Reconstruct pressure from the (merged) sons This must be specialised for each dime...
void further_build()
Further build for Crouzeix_Raviart interpolates the internal pressure dofs from father element: Make ...
double dshape_and_dtest_eulerian_nst(const Vector< double > &s, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const
Velocity shape and test functions and their derivs w.r.t. to global coords at local coordinate s (tak...
void pshape_nst(const Vector< double > &s, Shape &psi, Shape &test) const
Pressure shape and test functions at local coordinte s.
unsigned required_nvalue(const unsigned &n) const
Number of values (pinned or dofs) required at local node n.
void fix_pressure(const unsigned &p_dof, const double &p_value)
Pin p_dof-th pressure dof and set it to value specified by p_value.
PRefineableGeneralisedNewtonianQCrouzeixRaviartElement(const PRefineableGeneralisedNewtonianQCrouzeixRaviartElement< DIM > &dummy)=delete
Broken copy constructor.
unsigned ncont_interpolated_values() const
Number of continuously interpolated values: DIM (velocities)
void get_interpolated_values(const unsigned &t, const Vector< double > &s, Vector< double > &values)
Get all function values [u,v..,p] at previous timestep t (t=0: present; t>0: previous timestep).
void get_interpolated_values(const Vector< double > &s, Vector< double > &values)
Get the function value u in Vector. Note: Given the generality of the interface (this function is usu...
A class that is used to template the p-refineable Q elements by dimension. It's really nothing more t...
Definition: Qelements.h:2274
////////////////////////////////////////////////////////////////// //////////////////////////////////...
Definition: matrices.h:1370
RefineableElements are FiniteElements that may be subdivided into children to provide a better local ...
virtual RefineableElement * father_element_pt() const
Return a pointer to the father element.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
void dinterpolated_u_nst_ddata(const Vector< double > &s, const unsigned &i, Vector< double > &du_ddata, Vector< unsigned > &global_eqn_number)
Compute the derivatives of the i-th component of velocity at point s with respect to all data that ca...
virtual void get_dresidual_dnodal_coordinates(RankThreeTensor< double > &dresidual_dnodal_coordinates)
Compute derivatives of elemental residual vector with respect to nodal coordinates....
static void unpin_all_pressure_dofs(const Vector< GeneralisedElement * > &element_pt)
Unpin all pressure dofs in elements listed in vector.
static void pin_redundant_nodal_pressures(const Vector< GeneralisedElement * > &element_pt)
Loop over all elements in Vector (which typically contains all the elements in a fluid mesh) and pin ...
virtual void unpin_elemental_pressure_dofs()=0
Unpin all pressure dofs in the element.
virtual Node * pressure_node_pt(const unsigned &n_p)
Pointer to n_p-th pressure node (Default: NULL, indicating that pressure is not based on nodal interp...
void get_Z2_flux(const Vector< double > &s, Vector< double > &flux)
Get 'flux' for Z2 error recovery: Upper triangular entries in strain rate tensor.
virtual void pin_elemental_redundant_nodal_pressure_dofs()
Pin unused nodal pressure dofs (empty by default, because by default pressure dofs are not associated...
void fill_in_generic_residual_contribution_nst(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, unsigned flag)
Add element's contribution to elemental residual vector and/or Jacobian matrix flag=1: compute both f...
void get_pressure_and_velocity_mass_matrix_diagonal(Vector< double > &press_mass_diag, Vector< double > &veloc_mass_diag, const unsigned &which_one=0)
Compute the diagonal of the velocity/pressure mass matrices. If which one=0, both are computed,...
//////////////////////////////////////////////////////////////////////// ////////////////////////////...
RefineableGeneralisedNewtonianQCrouzeixRaviartElement(const RefineableGeneralisedNewtonianQCrouzeixRaviartElement< DIM > &dummy)=delete
Broken copy constructor.
void identify_load_data(std::set< std::pair< Data *, unsigned > > &paired_load_data)
Add to the set paired_load_data pairs containing.
void rebuild_from_sons(Mesh *&mesh_pt)
Rebuild from sons: Reconstruct pressure from the (merged) sons This must be specialised for each dime...
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation: Same order as shape functions.
void get_interpolated_values(const unsigned &t, const Vector< double > &s, Vector< double > &values)
Get all function values [u,v..,p] at previous timestep t (t=0: present; t>0: previous timestep).
void further_setup_hanging_nodes()
Perform additional hanging node procedures for variables that are not interpolated by all nodes....
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
void further_build()
Further build for Crouzeix_Raviart interpolates the internal pressure dofs from father element: Make ...
void get_interpolated_values(const Vector< double > &s, Vector< double > &values)
Get the function value u in Vector. Note: Given the generality of the interface (this function is usu...
Refineable version of Taylor Hood elements. These classes can be written in total generality.
unsigned ncont_interpolated_values() const
Number of continuously interpolated values: (DIM velocities + 1 pressure)
unsigned ninterpolating_node(const int &value_id)
The number of pressure nodes is 2^DIM. The number of velocity nodes is the same as the number of geom...
void identify_load_data(std::set< std::pair< Data *, unsigned > > &paired_load_data)
Add to the set paired_load_data pairs containing.
void further_setup_hanging_nodes()
Perform additional hanging node procedures for variables that are not interpolated by all nodes....
void get_interpolated_values(const unsigned &t, const Vector< double > &s, Vector< double > &values)
Get the function value u in Vector. Note: Given the generality of the interface (this function is usu...
Node * interpolating_node_pt(const unsigned &n, const int &value_id)
The velocities are isoparametric and so the "nodes" interpolating the velocities are the geometric no...
void pin_elemental_redundant_nodal_pressure_dofs()
Pin all nodal pressure dofs that are not required.
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 required_nvalue(const unsigned &n) const
Number of values required at local node n. In order to simplify matters, we allocate storage for pres...
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation: Same order as shape functions.
void get_interpolated_values(const Vector< double > &s, Vector< double > &values)
Get the function value u in Vector. Note: Given the generality of the interface (this function is usu...
unsigned ninterpolating_node_1d(const int &value_id)
The number of 1d pressure nodes is 2, the number of 1d velocity nodes is the same as the number of 1d...
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 ...
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 n_value==DIM, the fraction is the same as the 1d nod...
Refineable version of QElement<3,NNODE_1D>.
A class that is used to template the refineable Q elements by dimension. It's really nothing more tha...
Definition: Qelements.h:2259
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
A slight extension to the standard template vector class so that we can include "graceful" array rang...
Definition: Vector.h:58
//////////////////////////////////////////////////////////////////// ////////////////////////////////...