Time dependent part(s) of physics for element interiors.
72 &libMesh::FEMContext::interior_rate,
73 &libMesh::DiffContext::get_elem_solution_rate_derivative,
77 const unsigned int n_u_dofs = context.get_dof_indices(this->
_disp_vars.
u()).size();
79 const std::vector<libMesh::Real> &JxW =
80 this->
get_fe(context)->get_JxW();
83 libMesh::DenseSubVector<libMesh::Number> &Fu = context.get_elem_residual(this->
_disp_vars.
u());
84 libMesh::DenseSubVector<libMesh::Number> &Fv = context.get_elem_residual(this->
_disp_vars.
v());
85 libMesh::DenseSubVector<libMesh::Number> &Fw = context.get_elem_residual(this->
_disp_vars.
w());
87 libMesh::DenseSubMatrix<libMesh::Number>& Kuu = context.get_elem_jacobian(this->
_disp_vars.
u(),this->
_disp_vars.
u());
88 libMesh::DenseSubMatrix<libMesh::Number>& Kuv = context.get_elem_jacobian(this->
_disp_vars.
u(),this->
_disp_vars.
v());
89 libMesh::DenseSubMatrix<libMesh::Number>& Kuw = context.get_elem_jacobian(this->
_disp_vars.
u(),this->
_disp_vars.
w());
91 libMesh::DenseSubMatrix<libMesh::Number>& Kvu = context.get_elem_jacobian(this->
_disp_vars.
v(),this->
_disp_vars.
u());
92 libMesh::DenseSubMatrix<libMesh::Number>& Kvv = context.get_elem_jacobian(this->
_disp_vars.
v(),this->
_disp_vars.
v());
93 libMesh::DenseSubMatrix<libMesh::Number>& Kvw = context.get_elem_jacobian(this->
_disp_vars.
v(),this->
_disp_vars.
w());
95 libMesh::DenseSubMatrix<libMesh::Number>& Kwu = context.get_elem_jacobian(this->
_disp_vars.
w(),this->
_disp_vars.
u());
96 libMesh::DenseSubMatrix<libMesh::Number>& Kwv = context.get_elem_jacobian(this->
_disp_vars.
w(),this->
_disp_vars.
v());
97 libMesh::DenseSubMatrix<libMesh::Number>& Kww = context.get_elem_jacobian(this->
_disp_vars.
w(),this->
_disp_vars.
w());
99 unsigned int n_qpoints = context.get_element_qrule().n_points();
102 const std::vector<std::vector<libMesh::Real> >& dphi_dxi =
103 this->
get_fe(context)->get_dphidxi();
105 const std::vector<std::vector<libMesh::Real> >& dphi_deta =
106 this->
get_fe(context)->get_dphideta();
108 const libMesh::DenseSubVector<libMesh::Number>& u_coeffs = context.get_elem_solution( this->
_disp_vars.
u() );
109 const libMesh::DenseSubVector<libMesh::Number>& v_coeffs = context.get_elem_solution( this->
_disp_vars.
v() );
110 const libMesh::DenseSubVector<libMesh::Number>& w_coeffs = context.get_elem_solution( this->
_disp_vars.
w() );
112 const libMesh::DenseSubVector<libMesh::Number>& dudt_coeffs = context.get_elem_solution_rate( this->
_disp_vars.
u() );
113 const libMesh::DenseSubVector<libMesh::Number>& dvdt_coeffs = context.get_elem_solution_rate( this->
_disp_vars.
v() );
114 const libMesh::DenseSubVector<libMesh::Number>& dwdt_coeffs = context.get_elem_solution_rate( this->
_disp_vars.
w() );
117 const std::vector<libMesh::RealGradient>& dxdxi = this->
get_fe(context)->get_dxyzdxi();
118 const std::vector<libMesh::RealGradient>& dxdeta = this->
get_fe(context)->get_dxyzdeta();
120 const unsigned int dim = 2;
122 for (
unsigned int qp=0; qp != n_qpoints; qp++)
125 libMesh::Gradient grad_u, grad_v, grad_w;
126 libMesh::Gradient dgradu_dt, dgradv_dt, dgradw_dt;
128 for(
unsigned int d = 0; d < n_u_dofs; d++ )
130 libMesh::RealGradient u_gradphi( dphi_dxi[d][qp], dphi_deta[d][qp] );
131 grad_u += u_coeffs(d)*u_gradphi;
132 grad_v += v_coeffs(d)*u_gradphi;
133 grad_w += w_coeffs(d)*u_gradphi;
135 dgradu_dt += dudt_coeffs(d)*u_gradphi;
136 dgradv_dt += dvdt_coeffs(d)*u_gradphi;
137 dgradw_dt += dwdt_coeffs(d)*u_gradphi;
140 libMesh::RealGradient grad_x( dxdxi[qp](0), dxdeta[qp](0) );
141 libMesh::RealGradient grad_y( dxdxi[qp](1), dxdeta[qp](1) );
142 libMesh::RealGradient grad_z( dxdxi[qp](2), dxdeta[qp](2) );
146 libMesh::Real lambda_sq = 0;
149 grad_u, grad_v, grad_w,
150 a_cov, a_contra, A_cov, A_contra,
158 libMesh::Real jac = JxW[qp];
160 for (
unsigned int i=0; i != n_u_dofs; i++)
162 libMesh::RealGradient u_gradphi( dphi_dxi[i][qp], dphi_deta[i][qp] );
164 for(
unsigned int alpha = 0; alpha < dim; alpha++ )
166 for(
unsigned int beta = 0; beta < dim; beta++ )
169 const libMesh::Real factor = common_factor*tau(alpha,beta);
171 libMesh::Real u_diagterm = dgradu_dt(beta)*u_gradphi(alpha) + dgradu_dt(alpha)*u_gradphi(beta);
172 libMesh::Real v_diagterm = dgradv_dt(beta)*u_gradphi(alpha) + dgradv_dt(alpha)*u_gradphi(beta);
173 libMesh::Real w_diagterm = dgradw_dt(beta)*u_gradphi(alpha) + dgradw_dt(alpha)*u_gradphi(beta);
175 Fu(i) += factor*u_diagterm;
176 Fv(i) += factor*v_diagterm;
177 Fw(i) += factor*w_diagterm;
179 for(
unsigned int lambda = 0; lambda < dim; lambda++ )
181 for(
unsigned int mu = 0; mu < dim; mu++ )
183 const libMesh::Real C1 = common_factor*C(alpha,beta,lambda,mu);
185 const libMesh::Real gamma_u = 0.5*( dgradu_dt(lambda)*(grad_x(mu)+grad_u(mu)) +
186 (grad_x(lambda)+grad_u(lambda))*dgradu_dt(mu) );
188 const libMesh::Real gamma_v = 0.5*( dgradv_dt(lambda)*(grad_y(mu)+grad_v(mu)) +
189 (grad_y(lambda)+grad_v(lambda))*dgradv_dt(mu) );
191 const libMesh::Real gamma_w = 0.5*( dgradw_dt(lambda)*(grad_z(mu)+grad_w(mu)) +
192 (grad_z(lambda)+grad_w(lambda))*dgradw_dt(mu) );
194 const libMesh::Real x_term = C1*( (grad_x(beta)+grad_u(beta))*u_gradphi(alpha) +
195 (grad_x(alpha)+grad_u(alpha))*u_gradphi(beta) );
197 const libMesh::Real y_term = C1*( (grad_y(beta)+grad_v(beta))*u_gradphi(alpha) +
198 (grad_y(alpha)+grad_v(alpha))*u_gradphi(beta) );
200 const libMesh::Real z_term = C1*( (grad_z(beta)+grad_w(beta))*u_gradphi(alpha) +
201 (grad_z(alpha)+grad_w(alpha))*u_gradphi(beta) );
203 const libMesh::Real gamma_sum = gamma_u + gamma_v + gamma_w;
205 Fu(i) += x_term*gamma_sum;
207 Fv(i) += y_term*gamma_sum;
209 Fw(i) += z_term*gamma_sum;
216 if( compute_jacobian )
218 for (
unsigned int i=0; i != n_u_dofs; i++)
220 libMesh::RealGradient u_gradphi_i( dphi_dxi[i][qp], dphi_deta[i][qp] );
222 for (
unsigned int j=0; j != n_u_dofs; j++)
224 libMesh::RealGradient u_gradphi_j( dphi_dxi[j][qp], dphi_deta[j][qp] );
226 for(
unsigned int alpha = 0; alpha < dim; alpha++ )
228 for(
unsigned int beta = 0; beta < dim; beta++ )
232 libMesh::Real u_diagterm = dgradu_dt(beta)*u_gradphi_i(alpha) + dgradu_dt(alpha)*u_gradphi_i(beta);
233 libMesh::Real v_diagterm = dgradv_dt(beta)*u_gradphi_i(alpha) + dgradv_dt(alpha)*u_gradphi_i(beta);
234 libMesh::Real w_diagterm = dgradw_dt(beta)*u_gradphi_i(alpha) + dgradw_dt(alpha)*u_gradphi_i(beta);
236 const libMesh::Real factor = common_factor*tau(alpha,beta);
238 const libMesh::Real ddiagterm = u_gradphi_j(beta)*u_gradphi_i(alpha) + u_gradphi_j(alpha)*u_gradphi_i(beta);
240 Kuu(i,j) += factor*ddiagterm*context.get_elem_solution_rate_derivative();
241 Kvv(i,j) += factor*ddiagterm*context.get_elem_solution_rate_derivative();
242 Kww(i,j) += factor*ddiagterm*context.get_elem_solution_rate_derivative();
244 for(
unsigned int lambda = 0; lambda < dim; lambda++ )
246 for(
unsigned int mu = 0; mu < dim; mu++ )
248 const libMesh::Real dgamma_du = 0.5*( u_gradphi_j(lambda)*(grad_x(mu)+grad_u(mu)) +
249 (grad_x(lambda)+grad_u(lambda))*u_gradphi_j(mu) );
251 const libMesh::Real dgamma_dv = 0.5*( u_gradphi_j(lambda)*(grad_y(mu)+grad_v(mu)) +
252 (grad_y(lambda)+grad_v(lambda))*u_gradphi_j(mu) );
254 const libMesh::Real dgamma_dw = 0.5*( u_gradphi_j(lambda)*(grad_z(mu)+grad_w(mu)) +
255 (grad_z(lambda)+grad_w(lambda))*u_gradphi_j(mu) );
257 const libMesh::Real C1 = common_factor*C(alpha,beta,lambda,mu);
259 const libMesh::Real dfactor_du = C1*dgamma_du*context.get_elem_solution_derivative();
260 const libMesh::Real dfactor_dv = C1*dgamma_dv*context.get_elem_solution_derivative();
261 const libMesh::Real dfactor_dw = C1*dgamma_dw*context.get_elem_solution_derivative();
263 Kuu(i,j) += u_diagterm*dfactor_du;
264 Kuv(i,j) += u_diagterm*dfactor_dv;
265 Kuw(i,j) += u_diagterm*dfactor_dw;
267 Kvu(i,j) += v_diagterm*dfactor_du;
268 Kvv(i,j) += v_diagterm*dfactor_dv;
269 Kvw(i,j) += v_diagterm*dfactor_dw;
271 Kwu(i,j) += w_diagterm*dfactor_du;
272 Kwv(i,j) += w_diagterm*dfactor_dv;
273 Kww(i,j) += w_diagterm*dfactor_dw;
275 const libMesh::Real gamma_u = 0.5*( dgradu_dt(lambda)*(grad_x(mu)+grad_u(mu)) +
276 (grad_x(lambda)+grad_u(lambda))*dgradu_dt(mu) );
278 const libMesh::Real gamma_v = 0.5*( dgradv_dt(lambda)*(grad_y(mu)+grad_v(mu)) +
279 (grad_y(lambda)+grad_v(lambda))*dgradv_dt(mu) );
281 const libMesh::Real gamma_w = 0.5*( dgradw_dt(lambda)*(grad_z(mu)+grad_w(mu)) +
282 (grad_z(lambda)+grad_w(lambda))*dgradw_dt(mu) );
284 const libMesh::Real x_term = C1*( (grad_x(beta)+grad_u(beta))*u_gradphi_i(alpha) +
285 (grad_x(alpha)+grad_u(alpha))*u_gradphi_i(beta) );
287 const libMesh::Real y_term = C1*( (grad_y(beta)+grad_v(beta))*u_gradphi_i(alpha) +
288 (grad_y(alpha)+grad_v(alpha))*u_gradphi_i(beta) );
290 const libMesh::Real z_term = C1*( (grad_z(beta)+grad_w(beta))*u_gradphi_i(alpha) +
291 (grad_z(alpha)+grad_w(alpha))*u_gradphi_i(beta) );
293 const libMesh::Real gamma_sum = gamma_u + gamma_v + gamma_w;
297 const libMesh::Real dxterm_du = C1*( u_gradphi_j(beta)*u_gradphi_i(alpha)
298 + u_gradphi_j(alpha)*u_gradphi_i(beta) )*context.get_elem_solution_derivative();
300 const libMesh::Real dyterm_dv = dxterm_du;
301 const libMesh::Real dzterm_dw = dxterm_du;
303 Kuu(i,j) += gamma_sum*dxterm_du;
305 Kvv(i,j) += gamma_sum*dyterm_dv;
307 Kww(i,j) += gamma_sum*dzterm_dw;
309 const libMesh::Real dgamma_sum_du = 0.5*( u_gradphi_j(lambda)*(grad_x(mu)+grad_u(mu))*context.get_elem_solution_rate_derivative()
310 + dgradu_dt(lambda)*u_gradphi_j(mu)*context.get_elem_solution_derivative()
311 + u_gradphi_j(mu)*(grad_x(lambda)+grad_u(lambda))*context.get_elem_solution_rate_derivative()
312 + dgradu_dt(mu)*u_gradphi_j(lambda)*context.get_elem_solution_derivative() );
314 const libMesh::Real dgamma_sum_dv = 0.5*( u_gradphi_j(lambda)*(grad_y(mu)+grad_v(mu))*context.get_elem_solution_rate_derivative()
315 + dgradv_dt(lambda)*u_gradphi_j(mu)*context.get_elem_solution_derivative()
316 + u_gradphi_j(mu)*(grad_y(lambda)+grad_v(lambda))*context.get_elem_solution_rate_derivative()
317 + dgradv_dt(mu)*u_gradphi_j(lambda)*context.get_elem_solution_derivative() );
319 const libMesh::Real dgamma_sum_dw = 0.5*( u_gradphi_j(lambda)*(grad_z(mu)+grad_w(mu))*context.get_elem_solution_rate_derivative()
320 + dgradw_dt(lambda)*u_gradphi_j(mu)*context.get_elem_solution_derivative()
321 + u_gradphi_j(mu)*(grad_z(lambda)+grad_w(lambda))*context.get_elem_solution_rate_derivative()
322 + dgradw_dt(mu)*u_gradphi_j(lambda)*context.get_elem_solution_derivative() );
324 Kuu(i,j) += x_term*dgamma_sum_du;
325 Kuv(i,j) += x_term*dgamma_sum_dv;
326 Kuw(i,j) += x_term*dgamma_sum_dw;
328 Kvu(i,j) += y_term*dgamma_sum_du;
329 Kvv(i,j) += y_term*dgamma_sum_dv;
330 Kvw(i,j) += y_term*dgamma_sum_dw;
332 Kwu(i,j) += z_term*dgamma_sum_du;
333 Kwv(i,j) += z_term*dgamma_sum_dv;
334 Kww(i,j) += z_term*dgamma_sum_dw;
void mass_residual_impl(bool compute_jacobian, AssemblyContext &context, InteriorFuncType interior_solution, VarDerivType get_solution_deriv, libMesh::Real mu=1.0)
Implementation of mass_residual.
DisplacementFEVariables _disp_vars
void compute_metric_tensors(unsigned int qp, const libMesh::FEBase &elem, const AssemblyContext &context, const libMesh::Gradient &grad_u, const libMesh::Gradient &grad_v, const libMesh::Gradient &grad_w, libMesh::TensorValue< libMesh::Real > &a_cov, libMesh::TensorValue< libMesh::Real > &a_contra, libMesh::TensorValue< libMesh::Real > &A_cov, libMesh::TensorValue< libMesh::Real > &A_contra, libMesh::Real &lambda_sq)
void compute_stress_and_elasticity(unsigned int dim, const libMesh::TensorValue< libMesh::Real > &g_contra, const libMesh::TensorValue< libMesh::Real > &g_cov, const libMesh::TensorValue< libMesh::Real > &G_contra, const libMesh::TensorValue< libMesh::Real > &G_cov, libMesh::TensorValue< libMesh::Real > &stress, ElasticityTensor &C)
StressStrainLaw _stress_strain_law
libMesh::Real _lambda_factor
libMesh::Real _h0
Membrane thickness.
const libMesh::FEGenericBase< libMesh::Real > * get_fe(const AssemblyContext &context)