29 #include "grins_config.h"
33 #include "libmesh/getpot.h"
34 #include "libmesh/quadrature.h"
35 #include "libmesh/boundary_info.h"
36 #include "libmesh/fem_system.h"
44 if( !input.have_variable(
"Physics/"+physics_name+
"/pressure") )
46 std::cerr <<
"Error: Must input pressure for "+physics_name+
"." << std::endl
47 <<
" Please set Physics/"+physics_name+
"/pressure." << std::endl;
64 const unsigned int n_u_dofs = context.get_dof_indices(
_disp_vars.
u_var()).size();
66 const std::vector<libMesh::Real> &JxW =
67 this->
get_fe(context)->get_JxW();
69 const std::vector<std::vector<libMesh::Real> >& u_phi =
70 this->
get_fe(context)->get_phi();
72 libMesh::DenseSubVector<libMesh::Number> &Fu = context.get_elem_residual(
_disp_vars.
u_var());
73 libMesh::DenseSubVector<libMesh::Number> &Fv = context.get_elem_residual(
_disp_vars.
v_var());
74 libMesh::DenseSubVector<libMesh::Number> &Fw = context.get_elem_residual(
_disp_vars.
w_var());
85 unsigned int n_qpoints = context.get_element_qrule().n_points();
88 const std::vector<std::vector<libMesh::Real> >& dphi_dxi =
89 this->
get_fe(context)->get_dphidxi();
91 const std::vector<std::vector<libMesh::Real> >& dphi_deta =
92 this->
get_fe(context)->get_dphideta();
94 const libMesh::DenseSubVector<libMesh::Number>& u_coeffs = context.get_elem_solution(
_disp_vars.
u_var() );
95 const libMesh::DenseSubVector<libMesh::Number>& v_coeffs = context.get_elem_solution(
_disp_vars.
v_var() );
96 const libMesh::DenseSubVector<libMesh::Number>& w_coeffs = context.get_elem_solution(
_disp_vars.
w_var() );
98 const std::vector<libMesh::RealGradient>& dxdxi = this->
get_fe(context)->get_dxyzdxi();
99 const std::vector<libMesh::RealGradient>& dxdeta = this->
get_fe(context)->get_dxyzdeta();
101 for (
unsigned int qp=0; qp != n_qpoints; qp++)
104 libMesh::Real sqrt_a = sqrt( dxdxi[qp]*dxdxi[qp]*dxdeta[qp]*dxdeta[qp]
105 - dxdxi[qp]*dxdeta[qp]*dxdeta[qp]*dxdxi[qp] );
108 libMesh::Gradient grad_u, grad_v, grad_w;
109 for(
unsigned int d = 0; d < n_u_dofs; d++ )
111 libMesh::RealGradient u_gradphi( dphi_dxi[d][qp], dphi_deta[d][qp] );
112 grad_u += u_coeffs(d)*u_gradphi;
113 grad_v += v_coeffs(d)*u_gradphi;
114 grad_w += w_coeffs(d)*u_gradphi;
117 libMesh::RealGradient dudxi( grad_u(0), grad_v(0), grad_w(0) );
118 libMesh::RealGradient dudeta( grad_u(1), grad_v(1), grad_w(1) );
120 libMesh::RealGradient A_1 = dxdxi[qp] + dudxi;
121 libMesh::RealGradient A_2 = dxdeta[qp] + dudeta;
123 libMesh::RealGradient A_3 = A_1.cross(A_2);
131 libMesh::RealGradient traction =
_pressure/sqrt_a*A_3;
133 libMesh::Real jac = JxW[qp];
135 for (
unsigned int i=0; i != n_u_dofs; i++)
137 Fu(i) += traction(0)*u_phi[i][qp]*jac;
139 Fv(i) += traction(1)*u_phi[i][qp]*jac;
141 Fw(i) += traction(2)*u_phi[i][qp]*jac;
143 if( compute_jacobian )
145 for (
unsigned int j=0; j != n_u_dofs; j++)
147 libMesh::RealGradient u_gradphi( dphi_dxi[j][qp], dphi_deta[j][qp] );
149 const libMesh::Real dt0_dv =
_pressure/sqrt_a*(u_gradphi(0)*A_2(2) - A_1(2)*u_gradphi(1));
150 const libMesh::Real dt0_dw =
_pressure/sqrt_a*(A_1(1)*u_gradphi(1) - u_gradphi(0)*A_2(1));
152 const libMesh::Real dt1_du =
_pressure/sqrt_a*(A_1(2)*u_gradphi(1) - u_gradphi(0)*A_2(2));
153 const libMesh::Real dt1_dw =
_pressure/sqrt_a*(u_gradphi(0)*A_2(0) - A_1(0)*u_gradphi(1));
155 const libMesh::Real dt2_du =
_pressure/sqrt_a*(u_gradphi(0)*A_2(1) - A_1(1)*u_gradphi(1));
156 const libMesh::Real dt2_dv =
_pressure/sqrt_a*(A_1(0)*u_gradphi(1) - u_gradphi(0)*A_2(0));
158 Kuv(i,j) += dt0_dv*u_phi[i][qp]*jac;
159 Kuw(i,j) += dt0_dw*u_phi[i][qp]*jac;
161 Kvu(i,j) += dt1_du*u_phi[i][qp]*jac;
162 Kvw(i,j) += dt1_dw*u_phi[i][qp]*jac;
164 Kwu(i,j) += dt2_du*u_phi[i][qp]*jac;
165 Kwv(i,j) += dt2_dv*u_phi[i][qp]*jac;
virtual void set_parameter(libMesh::Number ¶m_variable, const GetPot &input, const std::string ¶m_name, libMesh::Number param_default)
Each subclass can simultaneously read a parameter value from.
SolidMechanicsFEVariables _disp_vars
VariableIndex v_var() const
virtual void element_time_derivative(bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
Time dependent part(s) of physics for element interiors.
VariableIndex w_var() const
void reset_pressure(libMesh::Real pressure_in)
VariableIndex u_var() const
virtual ~ElasticMembraneConstantPressure()
const libMesh::FEGenericBase< libMesh::Real > * get_fe(const AssemblyContext &context)
ElasticMembraneConstantPressure()