32 #include "libmesh/getpot.h"
33 #include "libmesh/quadrature.h"
34 #include "libmesh/fem_system.h"
35 #include "libmesh/elem.h"
39 template<
typename StressStrainLaw>
45 _is_compressible(is_compressible)
48 template<
typename StressStrainLaw>
51 InteriorFuncType interior_solution,
52 VarDerivType get_solution_deriv,
55 const unsigned int n_u_dofs = context.get_dof_indices(this->_disp_vars.u()).size();
57 const std::vector<libMesh::Real> &JxW =
58 this->get_fe(context)->get_JxW();
60 const std::vector<std::vector<libMesh::Real> >& u_phi =
61 this->get_fe(context)->get_phi();
64 libMesh::DenseSubVector<libMesh::Number> &Fu = context.get_elem_residual(this->_disp_vars.u());
65 libMesh::DenseSubVector<libMesh::Number> &Fv = context.get_elem_residual(this->_disp_vars.v());
66 libMesh::DenseSubVector<libMesh::Number> &Fw = context.get_elem_residual(this->_disp_vars.w());
68 libMesh::DenseSubMatrix<libMesh::Number>& Kuu = context.get_elem_jacobian(this->_disp_vars.u(),this->_disp_vars.u());
69 libMesh::DenseSubMatrix<libMesh::Number>& Kvv = context.get_elem_jacobian(this->_disp_vars.v(),this->_disp_vars.v());
70 libMesh::DenseSubMatrix<libMesh::Number>& Kww = context.get_elem_jacobian(this->_disp_vars.w(),this->_disp_vars.w());
72 unsigned int n_qpoints = context.get_element_qrule().n_points();
74 for (
unsigned int qp=0; qp != n_qpoints; qp++)
76 libMesh::Real jac = JxW[qp];
78 libMesh::Real u_ddot, v_ddot, w_ddot;
79 (context.*interior_solution)( this->_disp_vars.u(), qp, u_ddot );
80 (context.*interior_solution)( this->_disp_vars.v(), qp, v_ddot );
81 (context.*interior_solution)( this->_disp_vars.w(), qp, w_ddot );
83 for (
unsigned int i=0; i != n_u_dofs; i++)
85 libMesh::Real value = this->_rho*this->_A*u_phi[i][qp]*jac*mu;
86 Fu(i) += value*u_ddot;
87 Fv(i) += value*v_ddot;
88 Fw(i) += value*w_ddot;
90 if( compute_jacobian )
92 for (
unsigned int j=0; j != n_u_dofs; j++)
94 libMesh::Real jac_term = mu*this->_rho*this->_A*u_phi[i][qp]*u_phi[j][qp]*jac;
95 jac_term *= (context.*get_solution_deriv)();
106 template<
typename StressStrainLaw>
108 const libMesh::FEBase& elem,
110 const libMesh::Gradient& grad_u,
111 const libMesh::Gradient& grad_v,
112 const libMesh::Gradient& grad_w,
117 libMesh::Real& lambda_sq )
119 const std::vector<libMesh::RealGradient>& dxdxi = elem.get_dxyzdxi();
121 const std::vector<libMesh::Real>& dxidx = elem.get_dxidx();
122 const std::vector<libMesh::Real>& dxidy = elem.get_dxidy();
123 const std::vector<libMesh::Real>& dxidz = elem.get_dxidz();
125 libMesh::RealGradient dxi( dxidx[qp], dxidy[qp], dxidz[qp] );
127 libMesh::RealGradient dudxi( grad_u(0), grad_v(0), grad_w(0) );
131 a_cov(0,0) = dxdxi[qp]*dxdxi[qp];
137 A_cov(0,0) = (dxdxi[qp] + dudxi)*(dxdxi[qp] + dudxi);
141 a_contra(0,0) = 1/a_cov(0,0);
147 A_contra(0,0) = 1/A_cov(0,0);
150 if( _is_compressible )
152 libmesh_not_implemented();
158 lambda_sq = a_cov(0,0)/A_cov(0,0);
162 A_cov(1,1) = lambda_sq;
163 A_cov(2,2) = lambda_sq;
164 A_contra(1,1) = 1.0/lambda_sq;
165 A_contra(2,2) = 1.0/lambda_sq;
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)
Helper functions for parsing material properties.
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.