35 #include "libmesh/quadrature.h"
43 _velocity_penalty_x_index(0),
44 _velocity_penalty_y_index(0),
45 _velocity_penalty_z_index(0),
46 _velocity_penalty_base_x_index(0),
47 _velocity_penalty_base_y_index(0),
48 _velocity_penalty_base_z_index(0)
62 context.get_element_fe(this->_flow_vars.u_var())->get_xyz();
63 context.get_element_fe(this->_flow_vars.u_var())->get_phi();
72 std::string section =
"Physics/"+this->_physics_name+
"/output_vars";
74 std::string vel_penalty =
"vel_penalty";
75 if (this->_physics_name ==
"VelocityPenalty2")
78 if (this->_physics_name ==
"VelocityPenalty3")
81 if( input.have_variable(section) )
83 unsigned int n_vars = input.vector_variable_size(section);
85 for(
unsigned int v = 0; v < n_vars; v++ )
87 std::string name = input(section,
"DIE!",v);
89 if( name == std::string(
"velocity_penalty") )
91 _velocity_penalty_x_index =
94 _velocity_penalty_y_index =
97 _velocity_penalty_z_index =
100 else if( name == std::string(
"velocity_penalty_base") )
102 _velocity_penalty_base_x_index =
105 _velocity_penalty_base_y_index =
108 _velocity_penalty_base_z_index =
113 std::cerr <<
"Error: Invalue output_vars value for "+this->_physics_name << std::endl
114 <<
" Found " << name << std::endl
115 <<
" Acceptable values are: velocity_penalty" << std::endl
116 <<
" velocity_penalty_base" << std::endl;
130 #ifdef GRINS_USE_GRVY_TIMERS
131 this->_timer->BeginTimer(
"VelocityPenalty::element_time_derivative");
135 const std::vector<libMesh::Real> &JxW =
136 context.get_element_fe(this->_flow_vars.u_var())->get_JxW();
139 const std::vector<std::vector<libMesh::Real> >& u_phi =
140 context.get_element_fe(this->_flow_vars.u_var())->get_phi();
142 const std::vector<libMesh::Point>& u_qpoint =
143 context.get_element_fe(this->_flow_vars.u_var())->get_xyz();
146 const unsigned int n_u_dofs = context.get_dof_indices(this->_flow_vars.u_var()).size();
149 libMesh::DenseSubMatrix<libMesh::Number> &Kuu = context.get_elem_jacobian(this->_flow_vars.u_var(), this->_flow_vars.u_var());
150 libMesh::DenseSubMatrix<libMesh::Number> &Kuv = context.get_elem_jacobian(this->_flow_vars.u_var(), this->_flow_vars.v_var());
151 libMesh::DenseSubMatrix<libMesh::Number> &Kvu = context.get_elem_jacobian(this->_flow_vars.v_var(), this->_flow_vars.u_var());
152 libMesh::DenseSubMatrix<libMesh::Number> &Kvv = context.get_elem_jacobian(this->_flow_vars.v_var(), this->_flow_vars.v_var());
154 libMesh::DenseSubMatrix<libMesh::Number>* Kwu = NULL;
155 libMesh::DenseSubMatrix<libMesh::Number>* Kwv = NULL;
156 libMesh::DenseSubMatrix<libMesh::Number>* Kww = NULL;
157 libMesh::DenseSubMatrix<libMesh::Number>* Kuw = NULL;
158 libMesh::DenseSubMatrix<libMesh::Number>* Kvw = NULL;
160 libMesh::DenseSubVector<libMesh::Number> &Fu = context.get_elem_residual(this->_flow_vars.u_var());
161 libMesh::DenseSubVector<libMesh::Number> &Fv = context.get_elem_residual(this->_flow_vars.v_var());
162 libMesh::DenseSubVector<libMesh::Number>* Fw = NULL;
164 if( this->_dim == 3 )
166 Kuw = &context.get_elem_jacobian(this->_flow_vars.u_var(), this->_flow_vars.w_var());
167 Kvw = &context.get_elem_jacobian(this->_flow_vars.v_var(), this->_flow_vars.w_var());
169 Kwu = &context.get_elem_jacobian(this->_flow_vars.w_var(), this->_flow_vars.u_var());
170 Kwv = &context.get_elem_jacobian(this->_flow_vars.w_var(), this->_flow_vars.v_var());
171 Kww = &context.get_elem_jacobian(this->_flow_vars.w_var(), this->_flow_vars.w_var());
172 Fw = &context.get_elem_residual(this->_flow_vars.w_var());
175 unsigned int n_qpoints = context.get_element_qrule().n_points();
177 for (
unsigned int qp=0; qp != n_qpoints; qp++)
180 libMesh::Number u, v;
181 u = context.interior_value(this->_flow_vars.u_var(), qp);
182 v = context.interior_value(this->_flow_vars.v_var(), qp);
184 libMesh::NumberVectorValue U(u,v);
186 U(2) = context.interior_value(this->_flow_vars.w_var(), qp);
188 libMesh::NumberVectorValue F;
189 libMesh::NumberTensorValue dFdU;
190 libMesh::NumberTensorValue* dFdU_ptr =
191 compute_jacobian ? &dFdU : NULL;
192 if (!this->compute_force(u_qpoint[qp], context.time, U, F, dFdU_ptr))
195 const libMesh::Real jac = JxW[qp];
197 for (
unsigned int i=0; i != n_u_dofs; i++)
199 const libMesh::Number jac_i = jac * u_phi[i][qp];
204 if( this->_dim == 3 )
206 (*Fw)(i) += F(2)*jac_i;
209 if( compute_jacobian )
211 for (
unsigned int j=0; j != n_u_dofs; j++)
213 const libMesh::Number jac_ij = context.get_elem_solution_derivative() * jac_i * u_phi[j][qp];
214 Kuu(i,j) += jac_ij * dFdU(0,0);
215 Kuv(i,j) += jac_ij * dFdU(0,1);
216 Kvu(i,j) += jac_ij * dFdU(1,0);
217 Kvv(i,j) += jac_ij * dFdU(1,1);
219 if( this->_dim == 3 )
221 (*Kuw)(i,j) += jac_ij * dFdU(0,2);
222 (*Kvw)(i,j) += jac_ij * dFdU(1,2);
224 (*Kwu)(i,j) += jac_ij * dFdU(2,0);
225 (*Kwv)(i,j) += jac_ij * dFdU(2,1);
226 (*Kww)(i,j) += jac_ij * dFdU(2,2);
234 #ifdef GRINS_USE_GRVY_TIMERS
235 this->_timer->EndTimer(
"VelocityPenalty::element_time_derivative");
244 const libMesh::Point& point,
245 libMesh::Real& value )
247 libMesh::DenseVector<libMesh::Number> output_vec(3);
249 if( quantity_index == this->_velocity_penalty_x_index )
251 (*this->normal_vector_function)(point, context.time, output_vec);
253 value = output_vec(0);
255 else if( quantity_index == this->_velocity_penalty_y_index )
257 (*this->normal_vector_function)(point, context.time, output_vec);
259 value = output_vec(1);
261 else if( quantity_index == this->_velocity_penalty_z_index )
263 (*this->normal_vector_function)(point, context.time, output_vec);
265 value = output_vec(2);
267 else if( quantity_index == this->_velocity_penalty_base_x_index )
269 (*this->base_velocity_function)(point, context.time, output_vec);
271 value = output_vec(0);
273 else if( quantity_index == this->_velocity_penalty_base_y_index )
275 (*this->base_velocity_function)(point, context.time, output_vec);
277 value = output_vec(1);
279 else if( quantity_index == this->_velocity_penalty_base_z_index )
281 (*this->base_velocity_function)(point, context.time, output_vec);
283 value = output_vec(2);
INSTANTIATE_INC_NS_SUBCLASS(VelocityPenalty)
unsigned int register_quantity(std::string name)
Register quantity to be postprocessed.
virtual void element_time_derivative(bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
Time dependent part(s) of physics for element interiors.
virtual void init_context(AssemblyContext &context)
Initialize context for added physics variables.
virtual void register_postprocessing_vars(const GetPot &input, PostProcessedQuantities< libMesh::Real > &postprocessing)
Register postprocessing variables for VelocityPenalty.
virtual void compute_postprocessed_quantity(unsigned int quantity_index, const AssemblyContext &context, const libMesh::Point &point, libMesh::Real &value)