35 #include "libmesh/quadrature.h"
43 _parsed_velocity_source_x_index(0),
44 _parsed_velocity_source_y_index(0),
45 _parsed_velocity_source_z_index(0)
59 context.get_element_fe(this->_flow_vars.u_var())->get_xyz();
60 context.get_element_fe(this->_flow_vars.u_var())->get_phi();
69 std::string section =
"Physics/"+this->_physics_name+
"/output_vars";
71 std::string vel_source =
"vel_source";
73 if( input.have_variable(section) )
75 unsigned int n_vars = input.vector_variable_size(section);
77 for(
unsigned int v = 0; v < n_vars; v++ )
79 std::string name = input(section,
"DIE!",v);
81 if( name == std::string(
"parsed_velocity_source") )
83 _parsed_velocity_source_x_index =
86 _parsed_velocity_source_y_index =
89 _parsed_velocity_source_z_index =
94 std::cerr <<
"Error: Invalue output_vars value for "+this->_physics_name << std::endl
95 <<
" Found " << name << std::endl
96 <<
" Acceptable values are: parsed_velocity_source" << std::endl;
110 #ifdef GRINS_USE_GRVY_TIMERS
111 this->_timer->BeginTimer(
"ParsedVelocitySource::element_time_derivative");
115 const std::vector<libMesh::Real> &JxW =
116 context.get_element_fe(this->_flow_vars.u_var())->get_JxW();
119 const std::vector<std::vector<libMesh::Real> >& u_phi =
120 context.get_element_fe(this->_flow_vars.u_var())->get_phi();
122 const std::vector<libMesh::Point>& u_qpoint =
123 context.get_element_fe(this->_flow_vars.u_var())->get_xyz();
126 const unsigned int n_u_dofs = context.get_dof_indices(this->_flow_vars.u_var()).size();
129 libMesh::DenseSubMatrix<libMesh::Number> &Kuu = context.get_elem_jacobian(this->_flow_vars.u_var(), this->_flow_vars.u_var());
130 libMesh::DenseSubMatrix<libMesh::Number> &Kuv = context.get_elem_jacobian(this->_flow_vars.u_var(), this->_flow_vars.v_var());
131 libMesh::DenseSubMatrix<libMesh::Number> &Kvu = context.get_elem_jacobian(this->_flow_vars.v_var(), this->_flow_vars.u_var());
132 libMesh::DenseSubMatrix<libMesh::Number> &Kvv = context.get_elem_jacobian(this->_flow_vars.v_var(), this->_flow_vars.v_var());
134 libMesh::DenseSubMatrix<libMesh::Number>* Kwu = NULL;
135 libMesh::DenseSubMatrix<libMesh::Number>* Kwv = NULL;
136 libMesh::DenseSubMatrix<libMesh::Number>* Kww = NULL;
137 libMesh::DenseSubMatrix<libMesh::Number>* Kuw = NULL;
138 libMesh::DenseSubMatrix<libMesh::Number>* Kvw = NULL;
140 libMesh::DenseSubVector<libMesh::Number> &Fu = context.get_elem_residual(this->_flow_vars.u_var());
141 libMesh::DenseSubVector<libMesh::Number> &Fv = context.get_elem_residual(this->_flow_vars.v_var());
142 libMesh::DenseSubVector<libMesh::Number>* Fw = NULL;
144 if( this->_dim == 3 )
146 Kuw = &context.get_elem_jacobian(this->_flow_vars.u_var(), this->_flow_vars.w_var());
147 Kvw = &context.get_elem_jacobian(this->_flow_vars.v_var(), this->_flow_vars.w_var());
149 Kwu = &context.get_elem_jacobian(this->_flow_vars.w_var(), this->_flow_vars.u_var());
150 Kwv = &context.get_elem_jacobian(this->_flow_vars.w_var(), this->_flow_vars.v_var());
151 Kww = &context.get_elem_jacobian(this->_flow_vars.w_var(), this->_flow_vars.w_var());
152 Fw = &context.get_elem_residual(this->_flow_vars.w_var());
155 unsigned int n_qpoints = context.get_element_qrule().n_points();
157 for (
unsigned int qp=0; qp != n_qpoints; qp++)
159 libMesh::NumberVectorValue F;
160 libMesh::NumberTensorValue dFdU;
161 libMesh::NumberTensorValue* dFdU_ptr =
162 compute_jacobian ? &dFdU : NULL;
163 if (!this->compute_force(u_qpoint[qp], context.time, context, F, dFdU_ptr))
166 const libMesh::Real jac = JxW[qp];
168 for (
unsigned int i=0; i != n_u_dofs; i++)
170 const libMesh::Number jac_i = jac * u_phi[i][qp];
175 if( this->_dim == 3 )
177 (*Fw)(i) += F(2)*jac_i;
180 if( compute_jacobian )
182 for (
unsigned int j=0; j != n_u_dofs; j++)
184 const libMesh::Number jac_ij = context.get_elem_solution_derivative() * jac_i * u_phi[j][qp];
185 Kuu(i,j) += jac_ij * dFdU(0,0);
186 Kuv(i,j) += jac_ij * dFdU(0,1);
187 Kvu(i,j) += jac_ij * dFdU(1,0);
188 Kvv(i,j) += jac_ij * dFdU(1,1);
190 if( this->_dim == 3 )
192 (*Kuw)(i,j) += jac_ij * dFdU(0,2);
193 (*Kvw)(i,j) += jac_ij * dFdU(1,2);
195 (*Kwu)(i,j) += jac_ij * dFdU(2,0);
196 (*Kwv)(i,j) += jac_ij * dFdU(2,1);
197 (*Kww)(i,j) += jac_ij * dFdU(2,2);
205 #ifdef GRINS_USE_GRVY_TIMERS
206 this->_timer->EndTimer(
"ParsedVelocitySource::element_time_derivative");
215 const libMesh::Point& point,
216 libMesh::Real& value )
218 libMesh::DenseVector<libMesh::Number> output_vec(3);
220 if( quantity_index == this->_parsed_velocity_source_x_index )
222 (*this->velocity_source_function)(context, point, context.time, output_vec);
224 value = output_vec(0);
226 else if( quantity_index == this->_parsed_velocity_source_y_index )
228 (*this->velocity_source_function)(context, point, context.time, output_vec);
230 value = output_vec(1);
232 else if( quantity_index == this->_parsed_velocity_source_z_index )
234 (*this->velocity_source_function)(context, point, context.time, output_vec);
236 value = output_vec(2);
virtual void register_postprocessing_vars(const GetPot &input, PostProcessedQuantities< libMesh::Real > &postprocessing)
Register postprocessing variables for ParsedVelocitySource.
unsigned int register_quantity(std::string name)
Register quantity to be postprocessed.
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 compute_postprocessed_quantity(unsigned int quantity_index, const AssemblyContext &context, const libMesh::Point &point, libMesh::Real &value)
INSTANTIATE_INC_NS_SUBCLASS(ParsedVelocitySource)