30 #include "grins_config.h"
37 #include "libmesh/fem_context.h"
38 #include "libmesh/quadrature.h"
48 _p_pinning(input,physics_name),
49 _pin_pressure( input(
"Physics/"+
PhysicsNaming::stokes()+
"/pin_pressure", false ) )
64 _p_pinning.check_pin_location(system.get_mesh());
72 #ifdef GRINS_USE_GRVY_TIMERS
73 this->_timer->BeginTimer(
"Stokes::element_time_derivative");
77 const unsigned int n_u_dofs = context.get_dof_indices(this->_flow_vars.u()).size();
78 const unsigned int n_p_dofs = context.get_dof_indices(this->_press_var.p()).size();
81 libmesh_assert (n_u_dofs == context.get_dof_indices(this->_flow_vars.v()).size());
83 libmesh_assert (n_u_dofs == context.get_dof_indices(this->_flow_vars.w()).size());
89 const std::vector<libMesh::Real> &JxW =
90 context.get_element_fe(this->_flow_vars.u())->get_JxW();
94 const std::vector<std::vector<libMesh::RealGradient> >& u_gradphi =
95 context.get_element_fe(this->_flow_vars.u())->get_dphi();
98 const std::vector<std::vector<libMesh::Real> >& p_phi =
99 context.get_element_fe(this->_press_var.p())->get_phi();
107 libMesh::DenseSubMatrix<libMesh::Number> &Kuu = context.get_elem_jacobian(this->_flow_vars.u(), this->_flow_vars.u());
108 libMesh::DenseSubMatrix<libMesh::Number> &Kvv = context.get_elem_jacobian(this->_flow_vars.v(), this->_flow_vars.v());
109 libMesh::DenseSubMatrix<libMesh::Number>* Kww = NULL;
111 libMesh::DenseSubMatrix<libMesh::Number> &Kup = context.get_elem_jacobian(this->_flow_vars.u(), this->_press_var.p());
112 libMesh::DenseSubMatrix<libMesh::Number> &Kvp = context.get_elem_jacobian(this->_flow_vars.v(), this->_press_var.p());
113 libMesh::DenseSubMatrix<libMesh::Number>* Kwp = NULL;
115 libMesh::DenseSubVector<libMesh::Number> &Fu = context.get_elem_residual(this->_flow_vars.u());
116 libMesh::DenseSubVector<libMesh::Number> &Fv = context.get_elem_residual(this->_flow_vars.v());
117 libMesh::DenseSubVector<libMesh::Number>* Fw = NULL;
119 if( this->_dim == 3 )
121 Kww = &context.get_elem_jacobian(this->_flow_vars.w(), this->_flow_vars.w());
122 Kwp = &context.get_elem_jacobian(this->_flow_vars.w(), this->_press_var.p());
123 Fw = &context.get_elem_residual(this->_flow_vars.w());
132 unsigned int n_qpoints = context.get_element_qrule().n_points();
134 for (
unsigned int qp=0; qp != n_qpoints; qp++)
137 libMesh::Number p, u, v, w;
138 p = context.interior_value(this->_press_var.p(), qp);
139 u = context.interior_value(this->_flow_vars.u(), qp);
140 v = context.interior_value(this->_flow_vars.v(), qp);
142 w = context.interior_value(this->_flow_vars.w(), qp);
144 libMesh::Gradient grad_u, grad_v, grad_w;
145 grad_u = context.interior_gradient(this->_flow_vars.u(), qp);
146 grad_v = context.interior_gradient(this->_flow_vars.v(), qp);
148 grad_w = context.interior_gradient(this->_flow_vars.w(), qp);
150 libMesh::NumberVectorValue Uvec (u,v);
155 libMesh::Real _mu_qp = this->_mu(context, qp);
160 for (
unsigned int i=0; i != n_u_dofs; i++)
163 ( p*u_gradphi[i][qp](0)
164 -_mu_qp*(u_gradphi[i][qp]*grad_u) );
167 ( p*u_gradphi[i][qp](1)
168 -_mu_qp*(u_gradphi[i][qp]*grad_v) );
171 (*Fw)(i) += JxW[qp] *
172 ( p*u_gradphi[i][qp](2)
173 -_mu_qp*(u_gradphi[i][qp]*grad_w) );
176 if (compute_jacobian)
178 for (
unsigned int j=0; j != n_u_dofs; j++)
185 Kuu(i,j) += JxW[qp] * context.get_elem_solution_derivative() *
186 (-_mu_qp*(u_gradphi[i][qp]*u_gradphi[j][qp]));
188 Kvv(i,j) += JxW[qp] * context.get_elem_solution_derivative() *
189 (-_mu_qp*(u_gradphi[i][qp]*u_gradphi[j][qp]));
193 (*Kww)(i,j) += JxW[qp] * context.get_elem_solution_derivative() *
194 (-_mu_qp*(u_gradphi[i][qp]*u_gradphi[j][qp]));
199 for (
unsigned int j=0; j != n_p_dofs; j++)
201 Kup(i,j) += context.get_elem_solution_derivative() * JxW[qp]*u_gradphi[i][qp](0)*p_phi[j][qp];
202 Kvp(i,j) += context.get_elem_solution_derivative() * JxW[qp]*u_gradphi[i][qp](1)*p_phi[j][qp];
204 (*Kwp)(i,j) += context.get_elem_solution_derivative() * JxW[qp]*u_gradphi[i][qp](2)*p_phi[j][qp];
212 #ifdef GRINS_USE_GRVY_TIMERS
213 this->_timer->EndTimer(
"Stokes::element_time_derivative");
224 #ifdef GRINS_USE_GRVY_TIMERS
225 this->_timer->BeginTimer(
"Stokes::element_constraint");
229 const unsigned int n_u_dofs = context.get_dof_indices(this->_flow_vars.u()).size();
230 const unsigned int n_p_dofs = context.get_dof_indices(this->_press_var.p()).size();
236 const std::vector<libMesh::Real> &JxW =
237 context.get_element_fe(this->_flow_vars.u())->get_JxW();
241 const std::vector<std::vector<libMesh::RealGradient> >& u_gradphi =
242 context.get_element_fe(this->_flow_vars.u())->get_dphi();
245 const std::vector<std::vector<libMesh::Real> >& p_phi =
246 context.get_element_fe(this->_press_var.p())->get_phi();
252 libMesh::DenseSubMatrix<libMesh::Number> &Kpu = context.get_elem_jacobian(this->_press_var.p(), this->_flow_vars.u());
253 libMesh::DenseSubMatrix<libMesh::Number> &Kpv = context.get_elem_jacobian(this->_press_var.p(), this->_flow_vars.v());
254 libMesh::DenseSubMatrix<libMesh::Number>* Kpw = NULL;
256 libMesh::DenseSubVector<libMesh::Number> &Fp = context.get_elem_residual(this->_press_var.p());
258 if( this->_dim == 3 )
260 Kpw = &context.get_elem_jacobian(this->_press_var.p(), this->_flow_vars.w());
264 unsigned int n_qpoints = context.get_element_qrule().n_points();
265 for (
unsigned int qp=0; qp != n_qpoints; qp++)
268 libMesh::Gradient grad_u, grad_v, grad_w;
269 grad_u = context.interior_gradient(this->_flow_vars.u(), qp);
270 grad_v = context.interior_gradient(this->_flow_vars.v(), qp);
272 grad_w = context.interior_gradient(this->_flow_vars.w(), qp);
276 for (
unsigned int i=0; i != n_p_dofs; i++)
278 Fp(i) += JxW[qp] * p_phi[i][qp] *
279 (grad_u(0) + grad_v(1));
281 Fp(i) += JxW[qp] * p_phi[i][qp] *
284 if (compute_jacobian)
286 for (
unsigned int j=0; j != n_u_dofs; j++)
288 Kpu(i,j) += context.get_elem_solution_derivative() * JxW[qp]*p_phi[i][qp]*u_gradphi[j][qp](0);
289 Kpv(i,j) += context.get_elem_solution_derivative() * JxW[qp]*p_phi[i][qp]*u_gradphi[j][qp](1);
291 (*Kpw)(i,j) += context.get_elem_solution_derivative() * JxW[qp]*p_phi[i][qp]*u_gradphi[j][qp](2);
303 _p_pinning.pin_value( context, compute_jacobian, this->_press_var.p() );
307 #ifdef GRINS_USE_GRVY_TIMERS
308 this->_timer->EndTimer(
"Stokes::element_constraint");
321 const std::vector<libMesh::Real> &JxW =
322 context.get_element_fe(this->_flow_vars.u())->get_JxW();
326 const std::vector<std::vector<libMesh::Real> >& u_phi =
327 context.get_element_fe(this->_flow_vars.u())->get_phi();
330 const unsigned int n_u_dofs = context.get_dof_indices(this->_flow_vars.u()).size();
333 libMesh::DenseSubVector<libMesh::Real> &F_u = context.get_elem_residual(this->_flow_vars.u());
334 libMesh::DenseSubVector<libMesh::Real> &F_v = context.get_elem_residual(this->_flow_vars.v());
335 libMesh::DenseSubVector<libMesh::Real>* F_w = NULL;
337 libMesh::DenseSubMatrix<libMesh::Real> &M_uu = context.get_elem_jacobian(this->_flow_vars.u(), this->_flow_vars.u());
338 libMesh::DenseSubMatrix<libMesh::Real> &M_vv = context.get_elem_jacobian(this->_flow_vars.v(), this->_flow_vars.v());
339 libMesh::DenseSubMatrix<libMesh::Real>* M_ww = NULL;
341 if( this->_dim == 3 )
343 F_w = &context.get_elem_residual(this->_flow_vars.w());
344 M_ww = &context.get_elem_jacobian(this->_flow_vars.w(), this->_flow_vars.w());
347 unsigned int n_qpoints = context.get_element_qrule().n_points();
349 for (
unsigned int qp = 0; qp != n_qpoints; ++qp)
356 libMesh::Real u_dot, v_dot, w_dot = 0.0;
357 context.interior_rate(this->_flow_vars.u(), qp, u_dot);
358 context.interior_rate(this->_flow_vars.v(), qp, v_dot);
360 if( this->_dim == 3 )
361 context.interior_rate(this->_flow_vars.w(), qp, w_dot);
363 for (
unsigned int i = 0; i != n_u_dofs; ++i)
365 F_u(i) -= JxW[qp]*this->_rho*u_dot*u_phi[i][qp];
366 F_v(i) -= JxW[qp]*this->_rho*v_dot*u_phi[i][qp];
368 if( this->_dim == 3 )
369 (*F_w)(i) -= JxW[qp]*this->_rho*w_dot*u_phi[i][qp];
371 if( compute_jacobian )
373 for (
unsigned int j=0; j != n_u_dofs; j++)
377 libMesh::Real value = context.get_elem_solution_derivative() * JxW[qp]*this->_rho*u_phi[i][qp]*u_phi[j][qp];
384 (*M_ww)(i,j) -= value;
GRINS::ICHandlingBase * _ic_handler
Physics class for Incompressible Navier-Stokes.
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 auxiliary_init(MultiphysicsSystem &system)
Any auxillary initialization a Physics class may need.
virtual void element_time_derivative(bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
Time dependent part(s) of physics for element interiors.
Interface with libMesh for solving Multiphysics problems.
virtual void element_constraint(bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
Constraint part(s) of physics for element interiors.