30 #include "grins_config.h"
36 #include "libmesh/fem_context.h"
37 #include "libmesh/quadrature.h"
47 _p_pinning(input,physics_name),
48 _pin_pressure( input(
"Physics/"+
stokes+
"/pin_pressure", false ) )
68 #ifdef GRINS_USE_GRVY_TIMERS
69 this->_timer->BeginTimer(
"Stokes::element_time_derivative");
73 const unsigned int n_u_dofs = context.get_dof_indices(this->_flow_vars.u_var()).size();
74 const unsigned int n_p_dofs = context.get_dof_indices(this->_flow_vars.p_var()).size();
77 libmesh_assert (n_u_dofs == context.get_dof_indices(this->_flow_vars.v_var()).size());
79 libmesh_assert (n_u_dofs == context.get_dof_indices(this->_flow_vars.w_var()).size());
85 const std::vector<libMesh::Real> &JxW =
86 context.get_element_fe(this->_flow_vars.u_var())->get_JxW();
90 const std::vector<std::vector<libMesh::RealGradient> >& u_gradphi =
91 context.get_element_fe(this->_flow_vars.u_var())->get_dphi();
94 const std::vector<std::vector<libMesh::Real> >& p_phi =
95 context.get_element_fe(this->_flow_vars.p_var())->get_phi();
103 libMesh::DenseSubMatrix<libMesh::Number> &Kuu = context.get_elem_jacobian(this->_flow_vars.u_var(), this->_flow_vars.u_var());
104 libMesh::DenseSubMatrix<libMesh::Number> &Kvv = context.get_elem_jacobian(this->_flow_vars.v_var(), this->_flow_vars.v_var());
105 libMesh::DenseSubMatrix<libMesh::Number>* Kww = NULL;
107 libMesh::DenseSubMatrix<libMesh::Number> &Kup = context.get_elem_jacobian(this->_flow_vars.u_var(), this->_flow_vars.p_var());
108 libMesh::DenseSubMatrix<libMesh::Number> &Kvp = context.get_elem_jacobian(this->_flow_vars.v_var(), this->_flow_vars.p_var());
109 libMesh::DenseSubMatrix<libMesh::Number>* Kwp = NULL;
111 libMesh::DenseSubVector<libMesh::Number> &Fu = context.get_elem_residual(this->_flow_vars.u_var());
112 libMesh::DenseSubVector<libMesh::Number> &Fv = context.get_elem_residual(this->_flow_vars.v_var());
113 libMesh::DenseSubVector<libMesh::Number>* Fw = NULL;
115 if( this->_dim == 3 )
117 Kww = &context.get_elem_jacobian(this->_flow_vars.w_var(), this->_flow_vars.w_var());
118 Kwp = &context.get_elem_jacobian(this->_flow_vars.w_var(), this->_flow_vars.p_var());
119 Fw = &context.get_elem_residual(this->_flow_vars.w_var());
128 unsigned int n_qpoints = context.get_element_qrule().n_points();
130 for (
unsigned int qp=0; qp != n_qpoints; qp++)
133 libMesh::Number p, u, v, w;
134 p = context.interior_value(this->_flow_vars.p_var(), qp);
135 u = context.interior_value(this->_flow_vars.u_var(), qp);
136 v = context.interior_value(this->_flow_vars.v_var(), qp);
138 w = context.interior_value(this->_flow_vars.w_var(), qp);
140 libMesh::Gradient grad_u, grad_v, grad_w;
141 grad_u = context.interior_gradient(this->_flow_vars.u_var(), qp);
142 grad_v = context.interior_gradient(this->_flow_vars.v_var(), qp);
144 grad_w = context.interior_gradient(this->_flow_vars.w_var(), qp);
146 libMesh::NumberVectorValue Uvec (u,v);
151 libMesh::Real _mu_qp = this->_mu(context, qp);
156 for (
unsigned int i=0; i != n_u_dofs; i++)
159 ( p*u_gradphi[i][qp](0)
160 -_mu_qp*(u_gradphi[i][qp]*grad_u) );
163 ( p*u_gradphi[i][qp](1)
164 -_mu_qp*(u_gradphi[i][qp]*grad_v) );
167 (*Fw)(i) += JxW[qp] *
168 ( p*u_gradphi[i][qp](2)
169 -_mu_qp*(u_gradphi[i][qp]*grad_w) );
172 if (compute_jacobian)
174 for (
unsigned int j=0; j != n_u_dofs; j++)
181 Kuu(i,j) += JxW[qp] * context.get_elem_solution_derivative() *
182 (-_mu_qp*(u_gradphi[i][qp]*u_gradphi[j][qp]));
184 Kvv(i,j) += JxW[qp] * context.get_elem_solution_derivative() *
185 (-_mu_qp*(u_gradphi[i][qp]*u_gradphi[j][qp]));
189 (*Kww)(i,j) += JxW[qp] * context.get_elem_solution_derivative() *
190 (-_mu_qp*(u_gradphi[i][qp]*u_gradphi[j][qp]));
195 for (
unsigned int j=0; j != n_p_dofs; j++)
197 Kup(i,j) += context.get_elem_solution_derivative() * JxW[qp]*u_gradphi[i][qp](0)*p_phi[j][qp];
198 Kvp(i,j) += context.get_elem_solution_derivative() * JxW[qp]*u_gradphi[i][qp](1)*p_phi[j][qp];
200 (*Kwp)(i,j) += context.get_elem_solution_derivative() * JxW[qp]*u_gradphi[i][qp](2)*p_phi[j][qp];
208 #ifdef GRINS_USE_GRVY_TIMERS
209 this->_timer->EndTimer(
"Stokes::element_time_derivative");
220 #ifdef GRINS_USE_GRVY_TIMERS
221 this->_timer->BeginTimer(
"Stokes::element_constraint");
225 const unsigned int n_u_dofs = context.get_dof_indices(this->_flow_vars.u_var()).size();
226 const unsigned int n_p_dofs = context.get_dof_indices(this->_flow_vars.p_var()).size();
232 const std::vector<libMesh::Real> &JxW =
233 context.get_element_fe(this->_flow_vars.u_var())->get_JxW();
237 const std::vector<std::vector<libMesh::RealGradient> >& u_gradphi =
238 context.get_element_fe(this->_flow_vars.u_var())->get_dphi();
241 const std::vector<std::vector<libMesh::Real> >& p_phi =
242 context.get_element_fe(this->_flow_vars.p_var())->get_phi();
248 libMesh::DenseSubMatrix<libMesh::Number> &Kpu = context.get_elem_jacobian(this->_flow_vars.p_var(), this->_flow_vars.u_var());
249 libMesh::DenseSubMatrix<libMesh::Number> &Kpv = context.get_elem_jacobian(this->_flow_vars.p_var(), this->_flow_vars.v_var());
250 libMesh::DenseSubMatrix<libMesh::Number>* Kpw = NULL;
252 libMesh::DenseSubVector<libMesh::Number> &Fp = context.get_elem_residual(this->_flow_vars.p_var());
254 if( this->_dim == 3 )
256 Kpw = &context.get_elem_jacobian(this->_flow_vars.p_var(), this->_flow_vars.w_var());
260 unsigned int n_qpoints = context.get_element_qrule().n_points();
261 for (
unsigned int qp=0; qp != n_qpoints; qp++)
264 libMesh::Gradient grad_u, grad_v, grad_w;
265 grad_u = context.interior_gradient(this->_flow_vars.u_var(), qp);
266 grad_v = context.interior_gradient(this->_flow_vars.v_var(), qp);
268 grad_w = context.interior_gradient(this->_flow_vars.w_var(), qp);
272 for (
unsigned int i=0; i != n_p_dofs; i++)
274 Fp(i) += JxW[qp] * p_phi[i][qp] *
275 (grad_u(0) + grad_v(1));
277 Fp(i) += JxW[qp] * p_phi[i][qp] *
280 if (compute_jacobian)
282 for (
unsigned int j=0; j != n_u_dofs; j++)
284 Kpu(i,j) += context.get_elem_solution_derivative() * JxW[qp]*p_phi[i][qp]*u_gradphi[j][qp](0);
285 Kpv(i,j) += context.get_elem_solution_derivative() * JxW[qp]*p_phi[i][qp]*u_gradphi[j][qp](1);
287 (*Kpw)(i,j) += context.get_elem_solution_derivative() * JxW[qp]*p_phi[i][qp]*u_gradphi[j][qp](2);
299 _p_pinning.pin_value( context, compute_jacobian, this->_flow_vars.p_var() );
303 #ifdef GRINS_USE_GRVY_TIMERS
304 this->_timer->EndTimer(
"Stokes::element_constraint");
317 const std::vector<libMesh::Real> &JxW =
318 context.get_element_fe(this->_flow_vars.u_var())->get_JxW();
322 const std::vector<std::vector<libMesh::Real> >& u_phi =
323 context.get_element_fe(this->_flow_vars.u_var())->get_phi();
326 const unsigned int n_u_dofs = context.get_dof_indices(this->_flow_vars.u_var()).size();
329 libMesh::DenseSubVector<libMesh::Real> &F_u = context.get_elem_residual(this->_flow_vars.u_var());
330 libMesh::DenseSubVector<libMesh::Real> &F_v = context.get_elem_residual(this->_flow_vars.v_var());
331 libMesh::DenseSubVector<libMesh::Real>* F_w = NULL;
333 libMesh::DenseSubMatrix<libMesh::Real> &M_uu = context.get_elem_jacobian(this->_flow_vars.u_var(), this->_flow_vars.u_var());
334 libMesh::DenseSubMatrix<libMesh::Real> &M_vv = context.get_elem_jacobian(this->_flow_vars.v_var(), this->_flow_vars.v_var());
335 libMesh::DenseSubMatrix<libMesh::Real>* M_ww = NULL;
337 if( this->_dim == 3 )
339 F_w = &context.get_elem_residual(this->_flow_vars.w_var());
340 M_ww = &context.get_elem_jacobian(this->_flow_vars.w_var(), this->_flow_vars.w_var());
343 unsigned int n_qpoints = context.get_element_qrule().n_points();
345 for (
unsigned int qp = 0; qp != n_qpoints; ++qp)
352 libMesh::Real u_dot, v_dot, w_dot = 0.0;
353 context.interior_rate(this->_flow_vars.u_var(), qp, u_dot);
354 context.interior_rate(this->_flow_vars.v_var(), qp, v_dot);
356 if( this->_dim == 3 )
357 context.interior_rate(this->_flow_vars.w_var(), qp, w_dot);
359 for (
unsigned int i = 0; i != n_u_dofs; ++i)
361 F_u(i) -= JxW[qp]*this->_rho*u_dot*u_phi[i][qp];
362 F_v(i) -= JxW[qp]*this->_rho*v_dot*u_phi[i][qp];
364 if( this->_dim == 3 )
365 (*F_w)(i) -= JxW[qp]*this->_rho*w_dot*u_phi[i][qp];
367 if( compute_jacobian )
369 for (
unsigned int j=0; j != n_u_dofs; j++)
373 libMesh::Real value = context.get_elem_solution_derivative() * JxW[qp]*this->_rho*u_phi[i][qp]*u_phi[j][qp];
380 (*M_ww)(i,j) -= value;
GRINS::ICHandlingBase * _ic_handler
Physics class for Incompressible Navier-Stokes.
GRINS::BCHandlingBase * _bc_handler
INSTANTIATE_INC_NS_SUBCLASS(Stokes)
Base class for reading and handling initial conditions for physics classes.
virtual void mass_residual(bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
Mass matrix part(s) for element interiors. All boundary terms lie within the time_derivative part...
virtual void element_time_derivative(bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
Time dependent part(s) of physics for element interiors.
Base class for reading and handling boundary conditions for physics classes.
virtual void element_constraint(bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
Constraint part(s) of physics for element interiors.