34 #include "libmesh/quadrature.h"
35 #include "libmesh/boundary_info.h"
59 context.get_element_fe(this->_flow_vars.u_var())->get_xyz();
60 context.get_element_fe(this->_flow_vars.u_var())->get_phi();
71 #ifdef GRINS_USE_GRVY_TIMERS
72 this->_timer->BeginTimer(
"VelocityDrag::element_time_derivative");
76 const std::vector<libMesh::Real> &JxW =
77 context.get_element_fe(this->_flow_vars.u_var())->get_JxW();
80 const std::vector<std::vector<libMesh::Real> >& u_phi =
81 context.get_element_fe(this->_flow_vars.u_var())->get_phi();
83 const std::vector<libMesh::Point>& u_qpoint =
84 context.get_element_fe(this->_flow_vars.u_var())->get_xyz();
87 const unsigned int n_u_dofs = context.get_dof_indices(this->_flow_vars.u_var()).size();
90 libMesh::DenseSubMatrix<libMesh::Number> &Kuu = context.get_elem_jacobian(this->_flow_vars.u_var(), this->_flow_vars.u_var());
91 libMesh::DenseSubMatrix<libMesh::Number> &Kuv = context.get_elem_jacobian(this->_flow_vars.u_var(), this->_flow_vars.v_var());
92 libMesh::DenseSubMatrix<libMesh::Number> &Kvu = context.get_elem_jacobian(this->_flow_vars.v_var(), this->_flow_vars.u_var());
93 libMesh::DenseSubMatrix<libMesh::Number> &Kvv = context.get_elem_jacobian(this->_flow_vars.v_var(), this->_flow_vars.v_var());
95 libMesh::DenseSubMatrix<libMesh::Number>* Kwu = NULL;
96 libMesh::DenseSubMatrix<libMesh::Number>* Kwv = NULL;
97 libMesh::DenseSubMatrix<libMesh::Number>* Kww = NULL;
98 libMesh::DenseSubMatrix<libMesh::Number>* Kuw = NULL;
99 libMesh::DenseSubMatrix<libMesh::Number>* Kvw = NULL;
101 libMesh::DenseSubVector<libMesh::Number> &Fu = context.get_elem_residual(this->_flow_vars.u_var());
102 libMesh::DenseSubVector<libMesh::Number> &Fv = context.get_elem_residual(this->_flow_vars.v_var());
103 libMesh::DenseSubVector<libMesh::Number>* Fw = NULL;
105 if( this->_dim == 3 )
107 Kuw = &context.get_elem_jacobian(this->_flow_vars.u_var(), this->_flow_vars.w_var());
108 Kvw = &context.get_elem_jacobian(this->_flow_vars.v_var(), this->_flow_vars.w_var());
110 Kwu = &context.get_elem_jacobian(this->_flow_vars.w_var(), this->_flow_vars.u_var());
111 Kwv = &context.get_elem_jacobian(this->_flow_vars.w_var(), this->_flow_vars.v_var());
112 Kww = &context.get_elem_jacobian(this->_flow_vars.w_var(), this->_flow_vars.w_var());
113 Fw = &context.get_elem_residual(this->_flow_vars.w_var());
116 unsigned int n_qpoints = context.get_element_qrule().n_points();
118 for (
unsigned int qp=0; qp != n_qpoints; qp++)
121 libMesh::Number u, v;
122 u = context.interior_value(this->_flow_vars.u_var(), qp);
123 v = context.interior_value(this->_flow_vars.v_var(), qp);
125 libMesh::NumberVectorValue U(u,v);
127 U(2) = context.interior_value(this->_flow_vars.w_var(), qp);
129 libMesh::NumberVectorValue F;
130 libMesh::NumberTensorValue dFdU;
131 libMesh::NumberTensorValue* dFdU_ptr =
132 compute_jacobian ? &dFdU : NULL;
133 if (!this->compute_force(u_qpoint[qp], context.time, U, F, dFdU_ptr))
136 libMesh::Real jac = JxW[qp];
138 for (
unsigned int i=0; i != n_u_dofs; i++)
140 const libMesh::Number jac_i = jac * u_phi[i][qp];
145 if( this->_dim == 3 )
146 (*Fw)(i) += F(2)*jac_i;
148 if( compute_jacobian )
150 for (
unsigned int j=0; j != n_u_dofs; j++)
152 const libMesh::Number jac_ij =
153 jac_i * context.get_elem_solution_derivative() *
155 Kuu(i,j) += jac_ij * dFdU(0,0);
156 Kuv(i,j) += jac_ij * dFdU(0,1);
157 Kvu(i,j) += jac_ij * dFdU(1,0);
158 Kvv(i,j) += jac_ij * dFdU(1,1);
160 if( this->_dim == 3 )
162 (*Kuw)(i,j) += jac_ij * dFdU(0,2);
163 (*Kvw)(i,j) += jac_ij * dFdU(1,2);
164 (*Kwu)(i,j) += jac_ij * dFdU(2,0);
165 (*Kwv)(i,j) += jac_ij * dFdU(2,1);
166 (*Kww)(i,j) += jac_ij * dFdU(2,2);
174 #ifdef GRINS_USE_GRVY_TIMERS
175 this->_timer->EndTimer(
"VelocityDrag::element_time_derivative");
virtual void init_context(AssemblyContext &context)
Initialize context for added physics variables.
virtual void element_time_derivative(bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
Time dependent part(s) of physics for element interiors.
virtual void read_input_options(const GetPot &input)
Read options from GetPot input file.
INSTANTIATE_INC_NS_SUBCLASS(VelocityDrag)