32 #ifndef OOMPH_GENERALISED_NEWTONIAN_NAVIER_STOKES_ELEMENTS_HEADER 33 #define OOMPH_GENERALISED_NEWTONIAN_NAVIER_STOKES_ELEMENTS_HEADER 37 #include <oomph-lib-config.h> 42 #include "../generic/Qelements.h" 43 #include "../generic/fsi.h" 44 #include "../generic/projection.h" 45 #include "../generic/generalised_newtonian_constitutive_models.h" 96 const unsigned& which_one=0)=0;
133 template <
unsigned DIM>
143 typedef void (*NavierStokesBodyForceFctPt)(
const double& time,
149 typedef double (*NavierStokesSourceFctPt)(
const double& time,
156 typedef double (*NavierStokesPressureAdvDiffSourceFctPt)(
234 virtual double dshape_and_dtest_eulerian_at_knot_nst(
const unsigned &ipt,
243 virtual double dshape_and_dtest_eulerian_at_knot_nst(
256 virtual double dpshape_and_dptest_eulerian_nst(
const Vector<double> &
s,
260 DShape &dptestdx)
const=0;
274 if(Body_force_fct_pt == 0)
277 for(
unsigned i=0;
i<DIM;
i++) {result[
i] = 0.0;}
282 (*Body_force_fct_pt)(time,x,result);
304 get_body_force_nst(time,ipt,s,x,body_force);
310 for (
unsigned i=0;
i<DIM;
i++)
313 get_body_force_nst(time,ipt,s,x_pls,body_force_pls);
314 for (
unsigned j=0;j<DIM;j++)
316 d_body_force_dx(j,
i)=(body_force_pls[j]-body_force[j])/eps_fd;
336 if (Source_fct_pt == 0) {
return 0;}
338 else {
return (*Source_fct_pt)(time,x);}
358 double source=get_source_nst(time,ipt,x);
362 double source_pls=0.0;
364 for (
unsigned i=0;
i<DIM;
i++)
367 source_pls=get_source_nst(time,ipt,x_pls);
368 gradient[
i]=(source_pls-source)/eps_fd;
383 virtual void fill_in_generic_residual_contribution_nst(
391 virtual void fill_in_generic_dresidual_contribution_nst(
392 double*
const ¶meter_pt,
409 Body_force_fct_pt(0), Source_fct_pt(0),
412 Use_extrapolated_strainrate_to_compute_second_invariant(false),
413 ALE_is_disabled(false)
419 G_pt = &Default_Gravity_vector;
421 Viscosity_Ratio_pt = &Default_Physical_Ratio_Value;
422 Density_Ratio_pt = &Default_Physical_Ratio_Value;
434 const double &
re()
const {
return *Re_pt;}
437 const double &
re_st()
const {
return *ReSt_pt;}
460 const double &
re_invfr()
const {
return *ReInvFr_pt;}
473 {
return Body_force_fct_pt;}
477 {
return Body_force_fct_pt;}
487 {
return Constitutive_eqn_pt;}
492 Use_extrapolated_strainrate_to_compute_second_invariant=
false;
498 Use_extrapolated_strainrate_to_compute_second_invariant=
true;
502 virtual unsigned npres_nst()
const=0;
510 Shape &test)
const=0;
517 double u_nst(
const unsigned &n,
const unsigned &
i)
const 522 double u_nst(
const unsigned &
t,
const unsigned &n,
523 const unsigned &
i)
const 552 const unsigned u_nodal_index = this->u_index_nst(i);
555 const unsigned n_time = time_stepper_pt->
ntstorage();
557 for(
unsigned t=0;
t<n_time;
t++)
570 ALE_is_disabled=
true;
579 ALE_is_disabled=
false;
584 virtual double p_nst(
const unsigned &n_p)
const=0;
587 virtual double p_nst(
const unsigned &
t,
const unsigned &n_p)
const=0;
590 virtual void fix_pressure(
const unsigned &p_dof,
const double &p_value)=0;
598 double pressure_integral()
const;
602 void max_and_min_invariant_and_viscosity(
double& min_invariant,
603 double& max_invariant,
604 double& min_viscosity,
605 double& max_viscosity)
const;
608 double dissipation()
const;
617 double kin_energy()
const;
620 double d_kin_energy_dt()
const;
646 extrapolated_strain_rate(s,strain_rate);
680 get_traction(s,N,load);
689 const unsigned& which_one=0);
702 const unsigned& nplot)
const 710 for (
unsigned iplot=0;iplot<num_plot_points;iplot++)
717 if(i<DIM) {file_out << interpolated_u_nst(s,i) << std::endl;}
720 else if(i==DIM) {file_out << interpolated_p_nst(s) << std::endl;}
726 std::stringstream error_stream;
728 <<
"These Navier Stokes elements only store " << DIM+1 <<
" fields, " 729 <<
"but i is currently " << i << std::endl;
732 OOMPH_CURRENT_FUNCTION,
733 OOMPH_EXCEPTION_LOCATION);
757 std::stringstream error_stream;
759 <<
"These Navier Stokes elements only store " << DIM+1 <<
" fields,\n" 760 <<
"but i is currently " << i << std::endl;
763 OOMPH_CURRENT_FUNCTION,
764 OOMPH_EXCEPTION_LOCATION);
780 void output(std::ostream &outfile,
const unsigned &nplot);
792 void output(FILE* file_pt,
const unsigned &nplot);
800 full_output(outfile,nplot);
806 void full_output(std::ostream &outfile,
const unsigned &nplot);
812 void output_veloc(std::ostream &outfile,
const unsigned &nplot,
818 void output_vorticity(std::ostream &outfile,
819 const unsigned &nplot);
824 void output_fct(std::ostream &outfile,
const unsigned &nplot,
830 void output_fct(std::ostream &outfile,
const unsigned &nplot,
841 double& error,
double& norm);
849 double& error,
double& norm);
856 fill_in_generic_residual_contribution_nst(
867 fill_in_generic_residual_contribution_nst(residuals,jacobian,
880 fill_in_generic_residual_contribution_nst(residuals,jacobian,mass_matrix,2);
889 fill_in_generic_dresidual_contribution_nst(
902 fill_in_generic_dresidual_contribution_nst(
910 double*
const ¶meter_pt,
916 fill_in_generic_dresidual_contribution_nst(
917 parameter_pt,dres_dparam,djac_dparam,dmass_matrix_dparam,2);
928 for (
unsigned j=0;j<nint;j++)
931 if (eqn_number_backup[data_pt].
size()==0)
933 unsigned nvalue=data_pt->
nvalue();
934 eqn_number_backup[data_pt].resize(nvalue);
935 for (
unsigned i=0;
i<nvalue;
i++)
947 unsigned nnod=this->
nnode();
948 for (
unsigned j=0;j<nnod;j++)
952 if (eqn_number_backup[nod_pt].
size()==0)
955 unsigned nvalue=nod_pt->
nvalue();
956 eqn_number_backup[nod_pt].resize(nvalue);
957 for (
unsigned i=0;
i<nvalue;
i++)
987 if (eqn_number_backup[solid_posn_data_pt].
size()==0)
989 unsigned nvalue=solid_posn_data_pt->
nvalue();
990 eqn_number_backup[solid_posn_data_pt].resize(nvalue);
991 for (
unsigned i=0;
i<nvalue;
i++)
994 eqn_number_backup[solid_posn_data_pt][
i]=
998 solid_posn_data_pt->
pin(
i);
1012 dresidual_dnodal_coordinates);
1019 unsigned n_node =
nnode();
1025 for (
unsigned i=0;
i<DIM;
i++)
1028 unsigned u_nodal_index = u_index_nst(
i);
1032 for(
unsigned l=0;l<n_node;l++)
1043 unsigned n_node =
nnode();
1050 unsigned u_nodal_index = u_index_nst(i);
1053 double interpolated_u = 0.0;
1055 for(
unsigned l=0;l<n_node;l++)
1057 interpolated_u +=
nodal_value(l,u_nodal_index)*psi[l];
1060 return(interpolated_u);
1067 const unsigned &
i)
const 1070 unsigned n_node =
nnode();
1079 unsigned u_nodal_index = u_index_nst(i);
1082 double interpolated_u = 0.0;
1084 for(
unsigned l=0;l<n_node;l++)
1086 interpolated_u +=
nodal_value(t,l,u_nodal_index)*psi[l];
1089 return(interpolated_u);
1103 unsigned n_node =
nnode();
1110 const unsigned u_nodal_index = u_index_nst(i);
1114 for(
unsigned l=0;l<n_node;l++)
1118 if(global_eqn >= 0) {++n_u_dof;}
1122 du_ddata.resize(n_u_dof,0.0);
1123 global_eqn_number.resize(n_u_dof,0);
1128 for(
unsigned l=0;l<n_node;l++)
1132 if (global_eqn >= 0)
1135 global_eqn_number[count] = global_eqn;
1137 du_ddata[count] = psi[l];
1149 unsigned n_pres = npres_nst();
1156 double interpolated_p = 0.0;
1158 for(
unsigned l=0;l<n_pres;l++)
1160 interpolated_p += p_nst(l)*psi[l];
1163 return(interpolated_p);
1171 unsigned n_pres = npres_nst();
1178 double interpolated_p = 0.0;
1180 for(
unsigned l=0;l<n_pres;l++)
1182 interpolated_p += p_nst(t,l)*psi[l];
1185 return(interpolated_p);
1194 unsigned dim=s.size();
1197 data.resize(2*dim+1);
1200 for (
unsigned i=0;
i<
dim;
i++)
1203 data[
i+
dim]=this->interpolated_u_nst(s,
i);
1205 data[2*
dim]=this->interpolated_p_nst(s);
1221 template <
unsigned DIM>
1229 static const unsigned Initial_Nvalue[];
1250 inline double dshape_and_dtest_eulerian_at_knot_nst(
const unsigned &ipt,
1259 inline double dshape_and_dtest_eulerian_at_knot_nst(
1260 const unsigned &ipt,
1272 inline double dpshape_and_dptest_eulerian_nst(
const Vector<double> &
s,
1309 double p_nst(
const unsigned &
t,
const unsigned &
i)
const 1335 void identify_load_data(
1336 std::set<std::pair<Data*,unsigned> > &paired_load_data);
1346 void identify_pressure_data(
1347 std::set<std::pair<Data*,unsigned> > &paired_pressure_data);
1355 void output(std::ostream &outfile,
const unsigned &nplot)
1364 void output(FILE* file_pt,
const unsigned &nplot)
1395 std::list<std::pair<unsigned long,unsigned> >& dof_lookup_list)
const;
1406 template<
unsigned DIM>
1427 template<
unsigned DIM>
1430 const unsigned &ipt,
Shape &psi,
1469 for(
unsigned i=0;
i<9;
i++)
1473 for(
unsigned k=0;k<2;k++)
1475 dtestdx(
i,k) = dpsidx(
i,k);
1477 for(
unsigned p=0;p<2;p++)
1479 for(
unsigned q=0;q<9;q++)
1481 d_dtestdx_dX(p,q,
i,k) = d_dpsidx_dX(p,q,
i,k);
1517 for(
unsigned i=0;
i<27;
i++)
1521 for(
unsigned k=0;k<3;k++)
1523 dtestdx(
i,k) = dpsidx(
i,k);
1525 for(
unsigned p=0;p<3;p++)
1527 for(
unsigned q=0;q<27;q++)
1529 d_dtestdx_dX(p,q,
i,k) = d_dpsidx_dX(p,q,
i,k);
1616 template<
unsigned DIM>
1623 this->pshape_nst(s,psi);
1769 template <
unsigned DIM>
1776 static const unsigned Initial_Nvalue[];
1782 static const unsigned Pconv[];
1795 inline double dshape_and_dtest_eulerian_at_knot_nst(
const unsigned &ipt,
1804 inline double dshape_and_dtest_eulerian_at_knot_nst(
1805 const unsigned &ipt,
1818 inline double dpshape_and_dptest_eulerian_nst(
const Vector<double> &
s,
1833 {
return Initial_Nvalue[n];}
1857 double p_nst(
const unsigned &
t,
const unsigned &n_p)
const 1862 {
return static_cast<unsigned>(pow(2.0,static_cast<int>(DIM)));}
1879 void identify_load_data(
1880 std::set<std::pair<Data*,unsigned> > &paired_load_data);
1891 void identify_pressure_data(
1892 std::set<std::pair<Data*,unsigned> > &paired_pressure_data);
1900 void output(std::ostream &outfile,
const unsigned &nplot)
1908 void output(FILE* file_pt,
const unsigned &nplot)
1926 std::list<std::pair<unsigned long, unsigned> >& dof_lookup_list)
const;
1938 template<
unsigned DIM>
1963 template<
unsigned DIM>
1997 double psi1[2], psi2[2];
1998 double dpsi1[2],dpsi2[2];
2008 for(
unsigned i=0;
i<2;
i++)
2010 for(
unsigned j=0;j<2;j++)
2013 ppsi[2*
i+j] = psi2[
i]*psi1[j];
2014 dppsidx(2*
i+j,0)= psi2[
i]*dpsi1[j];
2015 dppsidx(2*
i+j,1)=dpsi2[
i]* psi1[j];
2071 for(
unsigned i=0;
i<9;
i++)
2075 for(
unsigned k=0;k<2;k++)
2077 dtestdx(
i,k) = dpsidx(
i,k);
2079 for(
unsigned p=0;p<2;p++)
2081 for(
unsigned q=0;q<9;q++)
2083 d_dtestdx_dX(p,q,
i,k) = d_dpsidx_dX(p,q,
i,k);
2119 for(
unsigned i=0;
i<27;
i++)
2123 for(
unsigned k=0;k<3;k++)
2125 dtestdx(
i,k) = dpsidx(
i,k);
2127 for(
unsigned p=0;p<3;p++)
2129 for(
unsigned q=0;q<27;q++)
2131 d_dtestdx_dX(p,q,
i,k) = d_dpsidx_dX(p,q,
i,k);
2155 double psi1[2], psi2[2];
2162 for(
unsigned i=0;
i<2;
i++)
2164 for(
unsigned j=0;j<2;j++)
2167 psi[2*
i + j] = psi2[
i]*psi1[j];
2188 double psi1[2], psi2[2], psi3[2];
2189 double dpsi1[2],dpsi2[2],dpsi3[2];
2201 for(
unsigned i=0;
i<2;
i++)
2203 for(
unsigned j=0;j<2;j++)
2205 for(
unsigned k=0;k<2;k++)
2208 ppsi[4*
i + 2*j + k] = psi3[
i]*psi2[j]*psi1[k];
2209 dppsidx(4*
i + 2*j + k,0) = psi3[
i] * psi2[j] *dpsi1[k];
2210 dppsidx(4*
i + 2*j + k,1) = psi3[
i] *dpsi2[j] * psi1[k];
2211 dppsidx(4*
i + 2*j + k,2) = dpsi3[
i] * psi2[j] * psi1[k];
2253 double psi1[2], psi2[2], psi3[2];
2262 for(
unsigned i=0;
i<2;
i++)
2264 for(
unsigned j=0;j<2;j++)
2266 for(
unsigned k=0;k<2;k++)
2269 psi[4*
i + 2*j + k] = psi3[
i]*psi2[j]*psi1[k];
2279 template<
unsigned DIM>
2287 this->pshape_nst(s,psi);
2354 template<
class TAYLOR_HOOD_ELEMENT>
2378 if (fld<this->
dim())
2381 unsigned nnod=this->
nnode();
2382 for (
unsigned j=0;j<nnod;j++)
2385 data_values.push_back(std::make_pair(this->
node_pt(j),fld));
2392 unsigned Pconv_size=this->
dim()+1;
2393 for (
unsigned j=0;j<Pconv_size;j++)
2396 unsigned vertex_index=this->Pconv[j];
2397 data_values.push_back(std::make_pair(this->
node_pt(vertex_index),fld));
2410 return this->
dim()+1;
2419 if (fld==this->
dim())
2443 unsigned n_dim=this->
dim();
2444 unsigned n_node=this->
nnode();
2449 this->pshape_nst(s,psi);
2451 Shape psif(n_node),testf(n_node);
2452 DShape dpsifdx(n_node,n_dim), dtestfdx(n_node,n_dim);
2455 double J=this->dshape_and_dtest_eulerian_nst(s,psif,dpsifdx,
2461 Shape testf(n_node);
2462 DShape dpsifdx(n_node,n_dim), dtestfdx(n_node,n_dim);
2465 double J=this->dshape_and_dtest_eulerian_nst(s,psi,dpsifdx,
2476 const unsigned &fld,
2479 unsigned n_dim =this->
dim();
2480 unsigned n_node=this->
nnode();
2485 return this->interpolated_p_nst(t,s);
2491 unsigned u_nodal_index = this->u_index_nst(fld);
2500 double interpolated_u = 0.0;
2503 for(
unsigned l=0;l<n_node;l++)
2505 interpolated_u += this->
nodal_value(t,l,u_nodal_index)*psi[l];
2507 return interpolated_u;
2516 if (fld==this->
dim())
2518 return this->npres_nst();
2522 return this->
nnode();
2531 if (fld==this->
dim())
2537 const unsigned u_nodal_index = this->u_index_nst(fld);
2549 template<
class ELEMENT>
2562 template<
class ELEMENT>
2575 template<
class CROUZEIX_RAVIART_ELEMENT>
2598 if (fld < this->
dim())
2601 const unsigned n_node=this->
nnode();
2602 for (
unsigned n=0;n<n_node;n++)
2605 data_values.push_back(std::make_pair(this->
node_pt(n),fld));
2612 const unsigned n_press = this->npres_nst();
2614 for(
unsigned j=0;j<n_press;j++)
2616 data_values.push_back(
2630 return this->
dim()+1;
2639 if (fld==this->
dim())
2663 unsigned n_dim=this->
dim();
2664 unsigned n_node=this->
nnode();
2669 this->pshape_nst(s,psi);
2671 Shape psif(n_node),testf(n_node);
2672 DShape dpsifdx(n_node,n_dim), dtestfdx(n_node,n_dim);
2675 double J=this->dshape_and_dtest_eulerian_nst(s,psif,dpsifdx,
2681 Shape testf(n_node);
2682 DShape dpsifdx(n_node,n_dim), dtestfdx(n_node,n_dim);
2685 double J=this->dshape_and_dtest_eulerian_nst(s,psi,dpsifdx,
2696 const unsigned &fld,
2699 unsigned n_dim =this->
dim();
2704 return this->interpolated_p_nst(s);
2709 return this->interpolated_u_nst(t,s,fld);
2717 if (fld==this->
dim())
2719 return this->npres_nst();
2723 return this->
nnode();
2732 if (fld==this->
dim())
2738 const unsigned u_nodal_index = this->u_index_nst(fld);
2750 template<
class ELEMENT>
2764 template<
class ELEMENT>
A Base class defining the generalise Newtonian constitutive relation.
int p_local_eqn(const unsigned &n) const
Return the local equation numbers for the pressure values.
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 ...
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...
virtual double get_source_nst(const double &time, const unsigned &ipt, const Vector< double > &x)
Calculate the source fct at given time and Eulerian position.
void get_load(const Vector< double > &s, const Vector< double > &N, Vector< double > &load)
This implements a pure virtual function defined in the FSIFluidElement class. The function computes t...
virtual double interpolated_p_nst(const Vector< double > &s) const
Return FE interpolated pressure at local coordinate s.
unsigned ndof_types() const
The number of "DOF types" that degrees of freedom in this element are sub-divided into: Velocity and ...
Vector< std::pair< Data *, unsigned > > data_values_of_field(const unsigned &fld)
Specify the values associated with field fld. The information is returned in a vector of pairs which ...
void output(std::ostream &outfile, const unsigned &nplot)
Redirect output to NavierStokesEquations output.
GeneralisedAxisymAdvectionDiffusionSourceFctPt Source_fct_pt
Pointer to source function:
double * Viscosity_Ratio_pt
Pointer to the viscosity ratio (relative to the viscosity used in the definition of the Reynolds numb...
double u_nst(const unsigned &t, const unsigned &n, const unsigned &i) const
Velocity i at local node n at timestep t (t=0: present; t>0: previous). Uses suitably interpolated va...
double jacobian_and_shape_of_field(const unsigned &fld, const Vector< double > &s, Shape &psi)
Return Jacobian of mapping and shape functions of field fld at local coordinate s.
virtual unsigned required_nvalue(const unsigned &n) const
Number of values that must be stored at local node n by the element. The default is 0...
void fill_in_contribution_to_jacobian_and_mass_matrix(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix)
Add the element's contribution to its residuals vector, jacobian matrix and mass matrix.
unsigned ndof_types() const
Returns the number of "DOF types" that degrees of freedom in this element are sub-divided into: Veloc...
double u_nst(const unsigned &n, const unsigned &i) const
Velocity i at local node n. Uses suitably interpolated value for hanging nodes. The use of u_index_ns...
virtual void get_s_plot(const unsigned &i, const unsigned &nplot, Vector< double > &s, const bool &shifted_to_interior=false) const
Get cector of local coordinates of plot point i (when plotting nplot points in each "coordinate direc...
virtual double dshape_eulerian_at_knot(const unsigned &ipt, Shape &psi, DShape &dpsidx) const
Return the geometric shape functions and also first derivatives w.r.t. global coordinates at the ipt-...
void dshape< 2 >(const double &s, double *DPsi)
Derivatives of 1D shape functions specialised to linear order (2 Nodes)
virtual void output(std::ostream &outfile)
Output the element data — typically the values at the nodes in a format suitable for post-processing...
double *& re_st_pt()
Pointer to product of Reynolds and Strouhal number (=Womersley number)
double dpshape_and_dptest_eulerian_nst(const Vector< double > &s, Shape &ppsi, DShape &dppsidx, Shape &ptest, DShape &dptestdx) const
Pressure shape and test functions and their derivs w.r.t. to global coords at local coordinate s (tak...
virtual void extrapolated_strain_rate(const unsigned &ipt, DenseMatrix< double > &strain_rate) const
Extrapolated strain-rate tensor: 1/2 (du_i/dx_j + du_j/dx_i) based on the previous time steps evaluat...
virtual unsigned nplot_points_paraview(const unsigned &nplot) const
Return the number of actual plot points for paraview plot with parameter nplot. Broken virtual; can b...
NavierStokesSourceFctPt & source_fct_pt()
Access function for the source-function pointer.
double p_nst(const unsigned &i) const
Return the i-th pressure value (Discontinous pressure interpolation – no need to cater for hanging 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.
virtual void fill_in_contribution_to_hessian_vector_products(Vector< double > const &Y, DenseMatrix< double > const &C, DenseMatrix< double > &product)
Fill in contribution to the product of the Hessian (derivative of Jacobian with respect to all variab...
static Vector< double > Gamma
Vector to decide whether the stress-divergence form is used or not.
void fill_in_contribution_to_djacobian_dparameter(double *const ¶meter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam)
Compute the element's residual Vector and the jacobian matrix Virtual function can be overloaded by h...
GeneralisedNewtonianConstitutiveEquation< DIM > * Constitutive_eqn_pt
Pointer to the generalised Newtonian constitutive equation.
double * ReSt_pt
Pointer to global Reynolds number x Strouhal number (=Womersley)
const double & viscosity_ratio() const
Viscosity ratio for element: Element's viscosity relative to the viscosity used in the definition of ...
void full_output(std::ostream &outfile)
Full output function: x,y,[z],u,v,[w],p,du/dt,dv/dt,[dw/dt],dissipation in tecplot format...
void use_extrapolated_strainrate_to_compute_second_invariant()
Use extrapolation for non-Newtonian const eqn.
void interpolated_u_nst(const Vector< double > &s, Vector< double > &veloc) const
Compute vector of FE interpolated velocity u at local coordinate s.
A general Finite Element class.
void output(std::ostream &outfile)
Redirect output to NavierStokesEquations output.
std::string scalar_name_paraview(const unsigned &i) const
Name of the i-th scalar field. Default implementation returns V1 for the first one, V2 for the second etc. Can (should!) be overloaded with more meaningful names in specific elements.
void(* UnsteadyExactSolutionFctPt)(const double &, const Vector< double > &, Vector< double > &)
Function pointer for function that computes Vector-valued time-dependent function as ...
TimeStepper *& position_time_stepper_pt()
Return a pointer to the position timestepper.
const double & density_ratio() const
Density ratio for element: Element's density relative to the viscosity used in the definition of the ...
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...
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 local coordinate s (tak...
const double & re_invfr() const
Global inverse Froude number.
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
void output(FILE *file_pt, const unsigned &nplot)
Redirect output to NavierStokesEquations output.
double * ReInvFr_pt
Pointer to global Reynolds number x inverse Froude number (= Bond number / Capillary number) ...
void output(FILE *file_pt, const unsigned &nplot)
Redirect output to NavierStokesEquations output.
virtual int p_nodal_index_nst() const
Set the value at which the pressure is stored in the nodes.
int local_equation(const unsigned &fld, const unsigned &j)
Return local equation number of value j in field fld.
unsigned P_nst_internal_index
long & eqn_number(const unsigned &i)
Return the equation number of the i-th stored variable.
unsigned nvalue() const
Return number of values stored in data object (incl pinned ones).
virtual unsigned u_index_nst(const unsigned &i) const
Return the index at which the i-th unknown velocity component.
void output(std::ostream &outfile)
Redirect output to NavierStokesEquations output.
static int Pressure_not_stored_at_node
Static "magic" number that indicates that the pressure is not stored at a node.
double interpolated_p_nst(const unsigned &t, const Vector< double > &s) const
Return FE interpolated pressure at local coordinate s at time level t.
unsigned nvalue_of_field(const unsigned &fld)
Return number of values in field fld.
void output(FILE *file_pt)
Redirect output to NavierStokesEquations output.
void pin(const unsigned &i)
Pin the i-th stored variable.
virtual 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_body_force_gradient_nst(const double &time, const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, DenseMatrix< double > &d_body_force_dx)
virtual void get_body_force_nst(const double &time, const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, Vector< double > &result)
Calculate the body force at a given time and local and/or Eulerian position. This function is virtual...
double du_dt_nst(const unsigned &n, const unsigned &i) const
i-th component of du/dt at local node n. Uses suitably interpolated value for hanging nodes...
unsigned nhistory_values_for_coordinate_projection()
Number of positional history values. (Note: count includes current value!)
double p_nst(const unsigned &n_p) const
Access function for the pressure values at local pressure node n_p (const version) ...
virtual double local_to_eulerian_mapping(const DShape &dpsids, DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
Calculate the mapping from local to Eulerian coordinates, given the derivatives of the shape function...
static bool Pre_multiply_by_viscosity_ratio
TimeStepper *& time_stepper_pt()
Access function for pointer to time stepper: Null if object is not time-dependent.
virtual ~GeneralisedNewtonianTemplateFreeNavierStokesEquationsBase()
Virtual destructor (empty)
double p_nst(const unsigned &t, const unsigned &n_p) const
Access function for the pressure values at local pressure node n_p (const version) ...
ProjectableGeneralisedNewtonianCrouzeixRaviartElement()
Constructor [this was only required explicitly from gcc 4.5.2 onwards...].
const Vector< double > & g() const
Vector of gravitational components.
void full_output(std::ostream &outfile, const unsigned &nplot)
Full output function: x,y,[z],u,v,[w],p,du/dt,dv/dt,[dw/dt],dissipation in tecplot format...
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as ...
unsigned npres_nst() const
Return number of pressure values.
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
void disable_ALE()
Disable ALE, i.e. assert the mesh is not moving – you do this at your own risk!
ProjectableGeneralisedNewtonianTaylorHoodElement()
Constructor [this was only required explicitly from gcc 4.5.2 onwards...].
GeneralisedNewtonianTemplateFreeNavierStokesEquationsBase()
Constructor (empty)
virtual int p_local_eqn(const unsigned &n) const =0
Access function for the local equation number information for the pressure. p_local_eqn[n] = local eq...
virtual void transform_derivatives(const DenseMatrix< double > &inverse_jacobian, DShape &dbasis) const
Convert derivative w.r.t.local coordinates to derivatives w.r.t the coordinates used to assemble the ...
double get_field(const unsigned &t, const unsigned &fld, const Vector< double > &s)
Return interpolated field fld at local coordinate s, at time level t (t=0: present; t>0: history valu...
int local_equation(const unsigned &fld, const unsigned &j)
Return local equation number of value j in field fld.
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 n...
GeneralisedNewtonianNavierStokesEquations()
Constructor: NULL the body force and source function and make sure the ALE terms are included by defa...
double dshape_eulerian(const Vector< double > &s, Shape &psi, DShape &dpsidx) const
Compute the geometric shape functions and also first derivatives w.r.t. global coordinates at local c...
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
unsigned nhistory_values_for_projection(const unsigned &fld)
Number of history values to be stored for fld-th field. Whatever the timestepper has set up for the v...
virtual void get_dof_numbers_for_unknowns(std::list< std::pair< unsigned long, unsigned > > &dof_lookup_list) const
Create a list of pairs for the unknowns that this element is "in charge of" – ignore any unknowns as...
unsigned n_u_nst() const
Return the number of velocity components Used in FluidInterfaceElements.
const double & re_st() const
Product of Reynolds and Strouhal number (=Womersley number)
Vector< double > * G_pt
Pointer to global gravity Vector.
void fill_in_contribution_to_djacobian_and_dmass_matrix_dparameter(double *const ¶meter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam, DenseMatrix< double > &dmass_matrix_dparam)
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
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...
Data *& internal_data_pt(const unsigned &i)
Return a pointer to i-th internal data object.
GeneralisedNewtonianQCrouzeixRaviartElement()
Constructor, there are DIM+1 internal values (for the pressure)
void point_output_data(const Vector< double > &s, Vector< double > &data)
Output solution in data vector at local cordinates s: x,y [,z], u,v,[w], p.
static Vector< double > Default_Gravity_vector
Static default value for the gravity vector.
Vector< std::pair< Data *, unsigned > > data_values_of_field(const unsigned &fld)
Specify the values associated with field fld. The information is returned in a vector of pairs which ...
NavierStokesSourceFctPt Source_fct_pt
Pointer to volumetric source function.
A class that represents a collection of data; each Data object may contain many different individual ...
virtual int p_nodal_index_nst() const
Return the index at which the pressure is stored if it is stored at the nodes. If not stored at the n...
virtual void dshape_local(const Vector< double > &s, Shape &psi, DShape &dpsids) const
Function to compute the geometric shape functions and derivatives w.r.t. local coordinates at local c...
unsigned npres_nst() const
Return number of pressure values.
Data *const & variable_position_pt() const
Pointer to variable_position data (const version)
virtual void output_fct(std::ostream &outfile, const unsigned &n_plot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output an exact solution over the element.
NavierStokesBodyForceFctPt & body_force_fct_pt()
Access function for the body-force pointer.
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
unsigned nscalar_paraview() const
Number of scalars/fields output by this element. Reimplements broken virtual function in base class...
double Default_Physical_Constant_Value
Default value for physical constants.
Vector< double > *& g_pt()
Pointer to Vector of gravitational components.
void pshape_nst(const Vector< double > &s, Shape &psi) const
Pressure shape functions at local coordinate s.
void enable_ALE()
(Re-)enable ALE, i.e. take possible mesh motion into account when evaluating the time-derivative. Note: By default, ALE is enabled, at the expense of possibly creating unnecessary work in problems where the mesh is, in fact, stationary.
void pin_all_non_pressure_dofs(std::map< Data *, std::vector< int > > &eqn_number_backup)
Pin all non-pressure dofs and backup eqn numbers.
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...
unsigned ntstorage() const
Return the number of doubles required to represent history (one for steady)
double *& re_pt()
Pointer to Reynolds number.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Compute the element's residual Vector.
unsigned ntstorage() const
Return total number of doubles stored per value to record time history of each value (one for steady ...
void fill_in_contribution_to_dresiduals_dparameter(double *const ¶meter_pt, Vector< double > &dres_dparam)
Compute the element's residual Vector.
double dpshape_and_dptest_eulerian_nst(const Vector< double > &s, Shape &ppsi, DShape &dppsidx, Shape &ptest, DShape &dptestdx) const
Pressure shape and test functions and their derivs w.r.t. to global coords at local coordinate s (tak...
void output(FILE *file_pt)
C-style output function: x,y,[z],u,v,[w],p in tecplot format. Default number of plot points...
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
void output(std::ostream &outfile, const unsigned &nplot)
Redirect output to NavierStokesEquations output.
unsigned ninternal_data() const
Return the number of internal data objects.
A Class for nodes that deform elastically (i.e. position is an unknown in the problem). The idea is that the Eulerian positions are stored in a Data object and the Lagrangian coordinates are stored in addition. The pointer that addresses the Eulerian positions is set to the pointer to Value in the Data object. Hence, SolidNode uses knowledge of the internal structure of Data and must be a friend of the Data class. In order to allow a mesh to deform via an elastic-style equation in deforming-domain problems, the positions are stored separately from the values, so that elastic problems may be combined with any other type of problem.
void scalar_value_paraview(std::ofstream &file_out, const unsigned &i, const unsigned &nplot) const
Write values of the i-th scalar field at the plot points. Needs to be implemented for each new specif...
void output(FILE *file_pt)
Redirect output to NavierStokesEquations output.
unsigned nfields_for_projection()
Number of fields to be projected: dim+1, corresponding to velocity components and pressure...
virtual unsigned required_nvalue(const unsigned &n) const
Number of values (pinned or dofs) required at node n. Can be overwritten for hanging node version...
void use_current_strainrate_to_compute_second_invariant()
Switch to fully implict evaluation of non-Newtonian const eqn.
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
bool ALE_is_disabled
Boolean flag to indicate if ALE formulation is disabled when time-derivatives are computed...
double interpolated_u_nst(const unsigned &t, const Vector< double > &s, const unsigned &i) const
Return FE interpolated velocity u[i] at local coordinate s at time level t (t=0: present; t>0: histor...
void pshape_nst(const Vector< double > &s, Shape &psi) const
Pressure shape functions at local coordinate s.
unsigned nvalue_of_field(const unsigned &fld)
Return number of values in field fld.
double get_field(const unsigned &t, const unsigned &fld, const Vector< double > &s)
Return interpolated field fld at local coordinate s, at time level t (t=0: present; t>0: history valu...
unsigned nhistory_values_for_coordinate_projection()
Number of positional history values (Note: count includes current value!)
unsigned add_internal_data(Data *const &data_pt, const bool &fd=true)
Add a (pointer to an) internal data object to the element and return the index required to obtain it ...
virtual void get_pressure_and_velocity_mass_matrix_diagonal(Vector< double > &press_mass_diag, Vector< double > &veloc_mass_diag, const unsigned &which_one=0)=0
Compute the diagonal of the velocity/pressure mass matrices. If which one=0, both are computed...
const double & re() const
Reynolds number.
double * Re_pt
Pointer to global Reynolds number.
bool is_steady() const
Flag to indicate if a timestepper has been made steady (possibly temporarily to switch off time-depen...
double jacobian_and_shape_of_field(const unsigned &fld, const Vector< double > &s, Shape &psi)
Return Jacobian of mapping and shape functions of field fld at local coordinate s.
Crouzeix Raviart upgraded to become projectable.
double *& density_ratio_pt()
Pointer to Density ratio.
virtual int p_nodal_index_nst() const =0
Return the index at which the pressure is stored if it is stored at the nodes. If not stored at the n...
unsigned nfields_for_projection()
Number of fields to be projected: dim+1, corresponding to velocity components and pressure...
void full_output(std::ostream &outfile)
Full output function: x,y,[z],u,v,[w],p,du/dt,dv/dt,[dw/dt],dissipation in tecplot format...
bool Use_extrapolated_strainrate_to_compute_second_invariant
unsigned nhistory_values_for_projection(const unsigned &fld)
Number of history values to be stored for fld-th field. Whatever the timestepper has set up for the v...
double interpolated_u_nst(const Vector< double > &s, const unsigned &i) const
Return FE interpolated velocity u[i] at local coordinate s.
static double Default_Physical_Constant_Value
Static default value for the physical constants (all initialised to zero)
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.
GeneralisedNewtonianQTaylorHoodElement()
Constructor, no internal data points.
virtual void get_source_gradient_nst(const double &time, const unsigned &ipt, const Vector< double > &x, Vector< double > &gradient)
double *& viscosity_ratio_pt()
Pointer to Viscosity Ratio.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Compute the element's residual Vector and the jacobian matrix Virtual function can be overloaded by h...
int p_local_eqn(const unsigned &n) const
Return the local equation numbers for the pressure values.
GeneralisedNewtonianConstitutiveEquation< DIM > *& constitutive_eqn_pt()
Access function for the constitutive equation pointer.
unsigned nnode() const
Return the number of nodes.
Base class for time-stepping schemes. Timestepper provides an approximation of the temporal derivativ...
virtual void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Plot the error when compared against a given exact solution . Also calculates the norm of the error a...
NavierStokesBodyForceFctPt body_force_fct_pt() const
Access function for the body-force pointer. Const version.
virtual double weight(const unsigned &i, const unsigned &j) const
Access function for j-th weight for the i-th derivative.
Taylor Hood upgraded to become projectable.
NavierStokesSourceFctPt source_fct_pt() const
Access function for the source-function pointer. Const version.
static double Default_fd_jacobian_step
Double used for the default finite difference step in elemental jacobian calculations.
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...
The FSIFluidElement class is a base class for all fluid finite elements that apply a load (traction) ...
static DenseMatrix< double > Dummy_matrix
Empty dense matrix used as a dummy argument to combined residual and jacobian functions in the case w...
double *& re_invfr_pt()
Pointer to global inverse Froude number.
std::string to_string(T object, unsigned float_precision=8)
Conversion function that should work for anything with operator<< defined (at least all basic types)...
double * Density_Ratio_pt
Pointer to the density ratio (relative to the density used in the definition of the Reynolds number) ...
static double Default_Physical_Ratio_Value
Static default value for the physical ratios (all are initialised to one)
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...
virtual void get_dresidual_dnodal_coordinates(RankThreeTensor< double > &dresidual_dnodal_coordinates)
Compute derivatives of elemental residual vector with respect to nodal coordinates. Default implementation by FD can be overwritten for specific elements. dresidual_dnodal_coordinates(l,i,j) = d res(l) / dX_{ij}.
virtual void pin_all_non_pressure_dofs(std::map< Data *, std::vector< int > > &eqn_number_backup)=0
Pin all non-pressure dofs and backup eqn numbers of all Data.
void output(std::ostream &outfile)
Output function: x,y,[z],u,v,[w],p in tecplot format. Default number of plot points.
NavierStokesBodyForceFctPt Body_force_fct_pt
Pointer to body force function.
int nodal_local_eqn(const unsigned &n, const unsigned &i) const
Return the local equation number corresponding to the i-th value at the n-th local node...
int internal_local_eqn(const unsigned &i, const unsigned &j) const
Return the local equation number corresponding to the j-th value stored at the i-th internal data...
void shape< 2 >(const double &s, double *Psi)
1D shape functions specialised to linear order (2 Nodes)