30 #include "libmesh/getpot.h"
31 #include "libmesh/mesh.h"
32 #include "libmesh/fem_system.h"
38 (
const std::string & helper_name,
45 if (input.have_variable(
"Stabilization/tau_constant_vel"))
47 (_C, input,
"Stabilization/tau_constant_vel", _C );
50 (_C, input,
"Stabilization/tau_constant", _C );
52 if (input.have_variable(
"Stabilization/tau_factor_vel"))
54 (_tau_factor, input,
"Stabilization/tau_factor_vel", _tau_factor );
57 (_tau_factor, input,
"Stabilization/tau_factor", _tau_factor );
75 libMesh::Gradient& grad_u,
76 libMesh::Gradient& grad_v )
const
78 return libMesh::RealGradient( U*grad_u, U*grad_v );
82 libMesh::Gradient& grad_u,
83 libMesh::Gradient& grad_v,
84 libMesh::Gradient& grad_w )
const
86 return libMesh::RealGradient( U*grad_u, U*grad_v, U*grad_w );
90 libMesh::RealTensor& hess_v )
const
92 return libMesh::RealGradient( hess_u(0,0) + hess_u(1,1),
93 hess_v(0,0) + hess_v(1,1) );
97 libMesh::RealTensor& hess_v,
98 libMesh::RealTensor& hess_w )
const
100 return libMesh::RealGradient( hess_u(0,0) + hess_u(1,1) + hess_u(2,2),
101 hess_v(0,0) + hess_v(1,1) + hess_v(2,2),
102 hess_w(0,0) + hess_w(1,1) + hess_w(2,2) );
106 libMesh::RealTensor& hess_v )
const
108 return libMesh::RealGradient( hess_u(0,0) + hess_v(0,1),
109 hess_u(1,0) + hess_v(1,1) );
113 libMesh::RealTensor& hess_v,
114 libMesh::RealTensor& hess_w )
const
116 return libMesh::RealGradient( hess_u(0,0) + hess_v(0,1) + hess_w(0,2),
117 hess_u(1,0) + hess_v(1,1) + hess_w(1,2),
118 hess_u(2,0) + hess_v(2,1) + hess_w(2,2) );
122 libMesh::RealTensor& hess_v )
const
124 return libMesh::RealGradient( hess_u(0,0) + hess_v(1,0),
125 hess_u(0,1) + hess_v(1,1) );
129 libMesh::RealTensor& hess_v,
130 libMesh::RealTensor& hess_w)
const
132 return libMesh::RealGradient( hess_u(0,0) + hess_v(1,0) + hess_w(2,0),
133 hess_u(0,1) + hess_v(1,1) + hess_w(2,1),
134 hess_u(0,2) + hess_v(1,2) + hess_w(2,2) );
138 unsigned int qp )
const
140 libMesh::RealGradient grad_u, grad_v;
142 grad_u = context.fixed_interior_gradient(this->
_flow_vars.
u_var(), qp);
143 grad_v = context.fixed_interior_gradient(this->
_flow_vars.
v_var(), qp);
145 libMesh::Real divU = grad_u(0) + grad_v(1);
147 if( context.get_system().get_mesh().mesh_dimension() == 3 )
149 divU += (context.fixed_interior_gradient(this->
_flow_vars.
w_var(), qp))(2);
158 libMesh::Real &res_C,
159 libMesh::Tensor &d_res_C_dgradU
162 libMesh::RealGradient grad_u, grad_v;
164 grad_u = context.fixed_interior_gradient(this->_flow_vars.u_var(), qp);
165 grad_v = context.fixed_interior_gradient(this->_flow_vars.v_var(), qp);
167 libMesh::Real divU = grad_u(0) + grad_v(1);
169 d_res_C_dgradU(0,0) = 1;
170 d_res_C_dgradU(1,1) = 1;
172 if( context.get_system().get_mesh().mesh_dimension() == 3 )
174 divU += (context.fixed_interior_gradient(this->_flow_vars.w_var(), qp))(2);
175 d_res_C_dgradU(2,2) = 1;
182 unsigned int qp,
const libMesh::Real rho,
const libMesh::Real mu )
const
184 libMesh::RealGradient U( context.fixed_interior_value(this->_flow_vars.u_var(), qp),
185 context.fixed_interior_value(this->_flow_vars.v_var(), qp) );
186 if(context.get_system().get_mesh().mesh_dimension() == 3)
189 libMesh::RealGradient grad_p = context.fixed_interior_gradient(this->
_flow_vars.
p_var(), qp);
191 libMesh::RealGradient grad_u = context.fixed_interior_gradient(this->
_flow_vars.
u_var(), qp);
192 libMesh::RealGradient grad_v = context.fixed_interior_gradient(this->
_flow_vars.
v_var(), qp);
194 libMesh::RealTensor hess_u = context.fixed_interior_hessian(this->
_flow_vars.
u_var(), qp);
195 libMesh::RealTensor hess_v = context.fixed_interior_hessian(this->
_flow_vars.
v_var(), qp);
197 libMesh::RealGradient rhoUdotGradU;
198 libMesh::RealGradient divGradU;
200 if( context.get_system().get_mesh().mesh_dimension() < 3 )
202 rhoUdotGradU = rho*this->
UdotGradU( U, grad_u, grad_v );
203 divGradU = this->
div_GradU( hess_u, hess_v );
207 libMesh::RealGradient grad_w = context.fixed_interior_gradient(this->
_flow_vars.
w_var(), qp);
208 libMesh::RealTensor hess_w = context.fixed_interior_hessian(this->
_flow_vars.
w_var(), qp);
210 rhoUdotGradU = rho*this->
UdotGradU( U, grad_u, grad_v, grad_w );
212 divGradU = this->
div_GradU( hess_u, hess_v, hess_w );
215 return rhoUdotGradU + grad_p - mu*divGradU;
220 unsigned int qp,
const libMesh::Real rho,
const libMesh::Real mu,
221 libMesh::Gradient &res_M,
222 libMesh::Tensor &d_res_M_dgradp,
223 libMesh::Tensor &d_res_M_dU,
224 libMesh::Gradient &d_res_Muvw_dgraduvw,
225 libMesh::Tensor &d_res_Muvw_dhessuvw
228 libMesh::RealGradient U( context.fixed_interior_value(this->_flow_vars.u_var(), qp),
229 context.fixed_interior_value(this->_flow_vars.v_var(), qp) );
230 if(context.get_system().get_mesh().mesh_dimension() == 3)
231 U(2) = context.fixed_interior_value(this->_flow_vars.w_var(), qp);
233 libMesh::RealGradient grad_p = context.fixed_interior_gradient(this->_flow_vars.p_var(), qp);
235 libMesh::RealGradient grad_u = context.fixed_interior_gradient(this->_flow_vars.u_var(), qp);
236 libMesh::RealGradient grad_v = context.fixed_interior_gradient(this->_flow_vars.v_var(), qp);
238 libMesh::RealTensor hess_u = context.fixed_interior_hessian(this->_flow_vars.u_var(), qp);
239 libMesh::RealTensor hess_v = context.fixed_interior_hessian(this->_flow_vars.v_var(), qp);
241 libMesh::RealGradient rhoUdotGradU;
242 libMesh::RealGradient divGradU;
244 if( context.get_system().get_mesh().mesh_dimension() < 3 )
246 rhoUdotGradU = rho*this->UdotGradU( U, grad_u, grad_v );
247 divGradU = this->div_GradU( hess_u, hess_v );
251 libMesh::RealGradient grad_w = context.fixed_interior_gradient(this->_flow_vars.w_var(), qp);
252 libMesh::RealTensor hess_w = context.fixed_interior_hessian(this->_flow_vars.w_var(), qp);
254 rhoUdotGradU = rho*this->UdotGradU( U, grad_u, grad_v, grad_w );
256 divGradU = this->div_GradU( hess_u, hess_v, hess_w );
258 d_res_M_dU(0,2) = rho * grad_u(2);
259 d_res_M_dU(1,2) = rho * grad_v(2);
261 d_res_Muvw_dgraduvw(2) = rho * U(2);
262 d_res_M_dgradp(2,2) = 1;
263 d_res_Muvw_dhessuvw(2,2) = mu;
265 d_res_M_dU(2,0) = rho * grad_w(0);
266 d_res_M_dU(2,1) = rho * grad_w(1);
267 d_res_M_dU(2,2) = rho * grad_w(2);
270 res_M = rhoUdotGradU + grad_p - mu*divGradU;
272 d_res_M_dU(0,0) = rho * grad_u(0);
273 d_res_M_dU(0,1) = rho * grad_u(1);
274 d_res_M_dU(1,0) = rho * grad_v(0);
275 d_res_M_dU(1,1) = rho * grad_v(1);
277 d_res_Muvw_dgraduvw(0) = rho * U(0);
278 d_res_Muvw_dgraduvw(1) = rho * U(1);
280 d_res_M_dgradp(0,0) = 1;
281 d_res_M_dgradp(1,1) = 1;
283 d_res_Muvw_dhessuvw(0,0) = mu;
284 d_res_Muvw_dhessuvw(1,1) = mu;
290 libMesh::RealGradient u_dot;
293 if(context.get_system().get_mesh().mesh_dimension() == 3)
294 context.interior_rate(this->_flow_vars.w_var(), qp, u_dot(2));
303 const libMesh::Real rho,
304 libMesh::RealGradient &res_M,
305 libMesh::Real &d_res_Muvw_duvw
308 libMesh::RealGradient u_dot;
309 context.interior_rate(this->_flow_vars.u_var(), qp, u_dot(0));
310 context.interior_rate(this->_flow_vars.v_var(), qp, u_dot(1));
311 if(context.get_system().get_mesh().mesh_dimension() == 3)
312 context.interior_rate(this->_flow_vars.w_var(), qp, u_dot(2));
315 d_res_Muvw_duvw = rho;
libMesh::RealGradient div_divU_I(libMesh::RealTensor &hess_u, libMesh::RealTensor &hess_v) const
VariableIndex w_var() const
libMesh::RealGradient compute_res_momentum_steady(AssemblyContext &context, unsigned int qp, const libMesh::Real rho, const libMesh::Real mu) const
libMesh::RealGradient compute_res_momentum_transient(AssemblyContext &context, unsigned int qp, const libMesh::Real rho) const
libMesh::RealGradient div_GradU(libMesh::RealTensor &hess_u, libMesh::RealTensor &hess_v) const
virtual void init(libMesh::FEMSystem *system)
void init(libMesh::FEMSystem &system)
VariableIndex p_var() const
PrimitiveFlowVariables _flow_vars
libMesh::RealGradient UdotGradU(libMesh::Gradient &U, libMesh::Gradient &grad_u, libMesh::Gradient &grad_v) const
VariableIndex u_var() const
libMesh::RealGradient div_GradU_T(libMesh::RealTensor &hess_u, libMesh::RealTensor &hess_v) const
void compute_res_momentum_steady_and_derivs(AssemblyContext &context, unsigned int qp, const libMesh::Real rho, const libMesh::Real mu, libMesh::Gradient &res_M, libMesh::Tensor &d_res_M_dgradp, libMesh::Tensor &d_res_M_dU, libMesh::Gradient &d_res_Muvw_dgraduvw, libMesh::Tensor &d_res_Muvw_dhessuvw) const
void compute_res_momentum_transient_and_derivs(AssemblyContext &context, unsigned int qp, const libMesh::Real rho, libMesh::RealGradient &res_M, libMesh::Real &d_res_Muvw_duvw) const
VariableIndex v_var() const
void compute_res_continuity_and_derivs(AssemblyContext &context, unsigned int qp, libMesh::Real &res_C, libMesh::Tensor &d_res_C_dgradU) const
IncompressibleNavierStokesStabilizationHelper(const std::string &helper_name, const GetPot &input)
~IncompressibleNavierStokesStabilizationHelper()
libMesh::Real compute_res_continuity(AssemblyContext &context, unsigned int qp) const