36 #include "libmesh/quadrature.h"
37 #include "libmesh/boundary_info.h"
45 _base_velocity_x_index(0),
46 _base_velocity_y_index(0),
47 _base_velocity_z_index(0)
62 context.get_element_fe(this->_flow_vars.u())->get_xyz();
63 context.get_element_fe(this->_flow_vars.u())->get_phi();
72 std::string section =
"Physics/"+this->_physics_name+
"/output_vars";
74 if( input.have_variable(section) )
76 unsigned int n_vars = input.vector_variable_size(section);
78 for(
unsigned int v = 0; v < n_vars; v++ )
80 std::string name = input(section,
"DIE!",v);
82 if( name == std::string(
"base_velocity") )
90 std::cerr <<
"Error: Invalid output_vars value for "+this->_physics_name << std::endl
91 <<
" Found " << name << std::endl
92 <<
" Acceptable values are: base_velocity" << std::endl;
104 #ifdef GRINS_USE_GRVY_TIMERS
105 this->_timer->BeginTimer(
"AveragedFan::element_time_derivative");
109 const std::vector<libMesh::Real> &JxW =
110 context.get_element_fe(this->_flow_vars.u())->get_JxW();
113 const std::vector<std::vector<libMesh::Real> >& u_phi =
114 context.get_element_fe(this->_flow_vars.u())->get_phi();
116 const std::vector<libMesh::Point>& u_qpoint =
117 context.get_element_fe(this->_flow_vars.u())->get_xyz();
120 const unsigned int n_u_dofs = context.get_dof_indices(this->_flow_vars.u()).size();
123 libMesh::DenseSubMatrix<libMesh::Number> &Kuu = context.get_elem_jacobian(this->_flow_vars.u(), this->_flow_vars.u());
124 libMesh::DenseSubMatrix<libMesh::Number> &Kuv = context.get_elem_jacobian(this->_flow_vars.u(), this->_flow_vars.v());
125 libMesh::DenseSubMatrix<libMesh::Number> &Kvu = context.get_elem_jacobian(this->_flow_vars.v(), this->_flow_vars.u());
126 libMesh::DenseSubMatrix<libMesh::Number> &Kvv = context.get_elem_jacobian(this->_flow_vars.v(), this->_flow_vars.v());
128 libMesh::DenseSubMatrix<libMesh::Number>* Kwu = NULL;
129 libMesh::DenseSubMatrix<libMesh::Number>* Kwv = NULL;
130 libMesh::DenseSubMatrix<libMesh::Number>* Kww = NULL;
131 libMesh::DenseSubMatrix<libMesh::Number>* Kuw = NULL;
132 libMesh::DenseSubMatrix<libMesh::Number>* Kvw = NULL;
134 libMesh::DenseSubVector<libMesh::Number> &Fu = context.get_elem_residual(this->_flow_vars.u());
135 libMesh::DenseSubVector<libMesh::Number> &Fv = context.get_elem_residual(this->_flow_vars.v());
136 libMesh::DenseSubVector<libMesh::Number>* Fw = NULL;
138 if( this->_dim == 3 )
140 Kuw = &context.get_elem_jacobian(this->_flow_vars.u(), this->_flow_vars.w());
141 Kvw = &context.get_elem_jacobian(this->_flow_vars.v(), this->_flow_vars.w());
143 Kwu = &context.get_elem_jacobian(this->_flow_vars.w(), this->_flow_vars.u());
144 Kwv = &context.get_elem_jacobian(this->_flow_vars.w(), this->_flow_vars.v());
145 Kww = &context.get_elem_jacobian(this->_flow_vars.w(), this->_flow_vars.w());
146 Fw = &context.get_elem_residual(this->_flow_vars.w());
149 unsigned int n_qpoints = context.get_element_qrule().n_points();
151 for (
unsigned int qp=0; qp != n_qpoints; qp++)
154 libMesh::Number u, v;
155 u = context.interior_value(this->_flow_vars.u(), qp);
156 v = context.interior_value(this->_flow_vars.v(), qp);
158 libMesh::NumberVectorValue U(u,v);
160 U(2) = context.interior_value(this->_flow_vars.w(), qp);
162 libMesh::NumberVectorValue F;
163 libMesh::NumberTensorValue dFdU;
164 libMesh::NumberTensorValue* dFdU_ptr =
165 compute_jacobian ? &dFdU : NULL;
166 if (!this->compute_force(u_qpoint[qp], context.time, U, F, dFdU_ptr))
169 libMesh::Real jac = JxW[qp];
171 for (
unsigned int i=0; i != n_u_dofs; i++)
173 const libMesh::Number jac_i = jac * u_phi[i][qp];
178 if( this->_dim == 3 )
179 (*Fw)(i) += F(2)*jac_i;
181 if( compute_jacobian )
183 for (
unsigned int j=0; j != n_u_dofs; j++)
185 const libMesh::Number jac_ij =
186 jac_i * context.get_elem_solution_derivative() *
188 Kuu(i,j) += jac_ij * dFdU(0,0);
189 Kuv(i,j) += jac_ij * dFdU(0,1);
190 Kvu(i,j) += jac_ij * dFdU(1,0);
191 Kvv(i,j) += jac_ij * dFdU(1,1);
193 if( this->_dim == 3 )
195 (*Kuw)(i,j) += jac_ij * dFdU(0,2);
196 (*Kvw)(i,j) += jac_ij * dFdU(1,2);
197 (*Kwu)(i,j) += jac_ij * dFdU(2,0);
198 (*Kwv)(i,j) += jac_ij * dFdU(2,1);
199 (*Kww)(i,j) += jac_ij * dFdU(2,2);
207 #ifdef GRINS_USE_GRVY_TIMERS
208 this->_timer->EndTimer(
"AveragedFan::element_time_derivative");
217 const libMesh::Point& point,
218 libMesh::Real& value )
220 libMesh::DenseVector<libMesh::Number> output_vec(3);
222 if( quantity_index == this->_base_velocity_x_index )
224 this->base_velocity_function(point, context.time, output_vec);
225 value = output_vec(0);
227 else if( quantity_index == this->_base_velocity_y_index )
229 this->base_velocity_function(point, context.time, output_vec);
230 value = output_vec(1);
232 else if( quantity_index == this->_base_velocity_z_index )
234 this->base_velocity_function(point, context.time, output_vec);
235 value = output_vec(2);
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.
INSTANTIATE_INC_NS_SUBCLASS(AveragedFan)
virtual void compute_postprocessed_quantity(unsigned int quantity_index, const AssemblyContext &context, const libMesh::Point &point, libMesh::Real &value)
Compute value of postprocessed quantities at libMesh::Point.
virtual void register_postprocessing_vars(const GetPot &input, PostProcessedQuantities< libMesh::Real > &postprocessing)
Register postprocessing variables for visualization output.
virtual void init_context(AssemblyContext &context)
Initialize context for added physics variables.