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>
42 bool is_compressible )
45 _is_compressible(is_compressible),
49 "Physics/"+physics_name+
"/h0",
61 template<
typename StressStrainLaw>
71 _lambda_sq_var = system->add_variable(
"lambda_sq", GRINSEnums::FIRST, GRINSEnums::LAGRANGE);
75 template<
typename StressStrainLaw>
78 InteriorFuncType interior_solution,
79 VarDerivType get_solution_deriv,
82 const unsigned int n_u_dofs = context.get_dof_indices(_disp_vars.u()).size();
84 const std::vector<libMesh::Real> &JxW =
85 this->get_fe(context)->get_JxW();
87 const std::vector<std::vector<libMesh::Real> >& u_phi =
88 this->get_fe(context)->get_phi();
91 libMesh::DenseSubVector<libMesh::Number> &Fu = context.get_elem_residual(_disp_vars.u());
92 libMesh::DenseSubVector<libMesh::Number> &Fv = context.get_elem_residual(_disp_vars.v());
93 libMesh::DenseSubVector<libMesh::Number> &Fw = context.get_elem_residual(_disp_vars.w());
95 libMesh::DenseSubMatrix<libMesh::Number>& Kuu = context.get_elem_jacobian(_disp_vars.u(),_disp_vars.u());
96 libMesh::DenseSubMatrix<libMesh::Number>& Kvv = context.get_elem_jacobian(_disp_vars.v(),_disp_vars.v());
97 libMesh::DenseSubMatrix<libMesh::Number>& Kww = context.get_elem_jacobian(_disp_vars.w(),_disp_vars.w());
99 unsigned int n_qpoints = context.get_element_qrule().n_points();
101 for (
unsigned int qp=0; qp != n_qpoints; qp++)
103 libMesh::Real jac = JxW[qp];
105 libMesh::Real u_ddot, v_ddot, w_ddot;
106 (context.*interior_solution)( _disp_vars.u(), qp, u_ddot );
107 (context.*interior_solution)( _disp_vars.v(), qp, v_ddot );
108 (context.*interior_solution)( _disp_vars.w(), qp, w_ddot );
110 for (
unsigned int i=0; i != n_u_dofs; i++)
112 Fu(i) += mu*this->_rho*_h0*u_ddot*u_phi[i][qp]*jac;
113 Fv(i) += mu*this->_rho*_h0*v_ddot*u_phi[i][qp]*jac;
114 Fw(i) += mu*this->_rho*_h0*w_ddot*u_phi[i][qp]*jac;
116 if( compute_jacobian )
118 for (
unsigned int j=0; j != n_u_dofs; j++)
120 libMesh::Real jac_term = this->_rho*_h0*u_phi[i][qp]*u_phi[j][qp]*jac;
121 jac_term *= mu*(context.*get_solution_deriv)();
123 Kuu(i,j) += jac_term;
124 Kvv(i,j) += jac_term;
125 Kww(i,j) += jac_term;
132 template<
typename StressStrainLaw>
134 const libMesh::FEBase& elem,
136 const libMesh::Gradient& grad_u,
137 const libMesh::Gradient& grad_v,
138 const libMesh::Gradient& grad_w,
143 libMesh::Real& lambda_sq )
145 const std::vector<libMesh::RealGradient>& dxdxi = elem.get_dxyzdxi();
146 const std::vector<libMesh::RealGradient>& dxdeta = elem.get_dxyzdeta();
148 const std::vector<libMesh::Real>& dxidx = elem.get_dxidx();
149 const std::vector<libMesh::Real>& dxidy = elem.get_dxidy();
150 const std::vector<libMesh::Real>& dxidz = elem.get_dxidz();
152 const std::vector<libMesh::Real>& detadx = elem.get_detadx();
153 const std::vector<libMesh::Real>& detady = elem.get_detady();
154 const std::vector<libMesh::Real>& detadz = elem.get_detadz();
156 libMesh::RealGradient dxi( dxidx[qp], dxidy[qp], dxidz[qp] );
157 libMesh::RealGradient deta( detadx[qp], detady[qp], detadz[qp] );
159 libMesh::RealGradient dudxi( grad_u(0), grad_v(0), grad_w(0) );
160 libMesh::RealGradient dudeta( grad_u(1), grad_v(1), grad_w(1) );
164 a_cov(0,0) = dxdxi[qp]*dxdxi[qp];
165 a_cov(0,1) = dxdxi[qp]*dxdeta[qp];
166 a_cov(1,0) = dxdeta[qp]*dxdxi[qp];
167 a_cov(1,1) = dxdeta[qp]*dxdeta[qp];
169 libMesh::Real det_a = a_cov(0,0)*a_cov(1,1) - a_cov(0,1)*a_cov(1,0);
173 A_cov(0,0) = (dxdxi[qp] + dudxi)*(dxdxi[qp] + dudxi);
174 A_cov(0,1) = (dxdxi[qp] + dudxi)*(dxdeta[qp] + dudeta);
175 A_cov(1,0) = (dxdeta[qp] + dudeta)*(dxdxi[qp] + dudxi);
176 A_cov(1,1) = (dxdeta[qp] + dudeta)*(dxdeta[qp] + dudeta);
180 a_contra(0,0) = dxi*dxi;
181 a_contra(0,1) = dxi*deta;
182 a_contra(1,0) = deta*dxi;
183 a_contra(1,1) = deta*deta;
186 libMesh::Real det_A = A_cov(0,0)*A_cov(1,1) - A_cov(0,1)*A_cov(1,0);
189 A_contra(0,0) = A_cov(1,1)/det_A;
190 A_contra(0,1) = -A_cov(0,1)/det_A;
191 A_contra(1,0) = -A_cov(1,0)/det_A;
192 A_contra(1,1) = A_cov(0,0)/det_A;
199 if( _is_compressible )
201 lambda_sq = context.interior_value(this->_lambda_sq_var, qp);
206 lambda_sq = det_a/det_A;
209 A_cov(2,2) = lambda_sq;
210 A_contra(2,2) = 1.0/lambda_sq;
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.
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)
libMesh::Real _rho
Membrane density.
static void read_density(const std::string &core_physics_name, const GetPot &input, ParameterUser ¶ms, libMesh::Real &rho)
Helper function to reading density from input.
static PhysicsName elastic_membrane()
Helper functions for parsing material properties.
virtual void init_variables(libMesh::FEMSystem *system)
Initialize variables for this physics.
static void read_property(const GetPot &input, const std::string &old_option, const std::string &property, const std::string &core_physics, ParameterUser ¶m_user, libMesh::Real &value)
Helper function for parsing/maintaing backward compatibility.
libMesh::Real _h0
Membrane thickness.
virtual void init_variables(libMesh::FEMSystem *system)
Initialize variables for this physics.