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;
57 libmesh_error_msg(
"ERROR: Cannot specify subdomains for "
59 +
"! Must specify subdomains through "
69 const unsigned int n_u_dofs = context.get_dof_indices(
_disp_vars.
u()).size();
71 const std::vector<libMesh::Real> &JxW =
72 this->
get_fe(context)->get_JxW();
74 const std::vector<std::vector<libMesh::Real> >& u_phi =
75 this->
get_fe(context)->get_phi();
77 libMesh::DenseSubVector<libMesh::Number> &Fu = context.get_elem_residual(
_disp_vars.
u());
78 libMesh::DenseSubVector<libMesh::Number> &Fv = context.get_elem_residual(
_disp_vars.
v());
79 libMesh::DenseSubVector<libMesh::Number> &Fw = context.get_elem_residual(
_disp_vars.
w());
90 unsigned int n_qpoints = context.get_element_qrule().n_points();
93 const std::vector<std::vector<libMesh::Real> >& dphi_dxi =
94 this->
get_fe(context)->get_dphidxi();
96 const std::vector<std::vector<libMesh::Real> >& dphi_deta =
97 this->
get_fe(context)->get_dphideta();
99 const libMesh::DenseSubVector<libMesh::Number>& u_coeffs = context.get_elem_solution(
_disp_vars.
u() );
100 const libMesh::DenseSubVector<libMesh::Number>& v_coeffs = context.get_elem_solution(
_disp_vars.
v() );
101 const libMesh::DenseSubVector<libMesh::Number>& w_coeffs = context.get_elem_solution(
_disp_vars.
w() );
103 const std::vector<libMesh::RealGradient>& dxdxi = this->
get_fe(context)->get_dxyzdxi();
104 const std::vector<libMesh::RealGradient>& dxdeta = this->
get_fe(context)->get_dxyzdeta();
106 for (
unsigned int qp=0; qp != n_qpoints; qp++)
109 libMesh::Real sqrt_a = sqrt( dxdxi[qp]*dxdxi[qp]*dxdeta[qp]*dxdeta[qp]
110 - dxdxi[qp]*dxdeta[qp]*dxdeta[qp]*dxdxi[qp] );
113 libMesh::Gradient grad_u, grad_v, grad_w;
114 for(
unsigned int d = 0; d < n_u_dofs; d++ )
116 libMesh::RealGradient u_gradphi( dphi_dxi[d][qp], dphi_deta[d][qp] );
117 grad_u += u_coeffs(d)*u_gradphi;
118 grad_v += v_coeffs(d)*u_gradphi;
119 grad_w += w_coeffs(d)*u_gradphi;
122 libMesh::RealGradient dudxi( grad_u(0), grad_v(0), grad_w(0) );
123 libMesh::RealGradient dudeta( grad_u(1), grad_v(1), grad_w(1) );
125 libMesh::RealGradient A_1 = dxdxi[qp] + dudxi;
126 libMesh::RealGradient A_2 = dxdeta[qp] + dudeta;
128 libMesh::RealGradient A_3 = A_1.cross(A_2);
136 libMesh::RealGradient traction =
_pressure/sqrt_a*A_3;
138 libMesh::Real jac = JxW[qp];
140 for (
unsigned int i=0; i != n_u_dofs; i++)
142 Fu(i) -= traction(0)*u_phi[i][qp]*jac;
144 Fv(i) -= traction(1)*u_phi[i][qp]*jac;
146 Fw(i) -= traction(2)*u_phi[i][qp]*jac;
148 if( compute_jacobian )
150 for (
unsigned int j=0; j != n_u_dofs; j++)
152 libMesh::RealGradient u_gradphi( dphi_dxi[j][qp], dphi_deta[j][qp] );
154 const libMesh::Real dt0_dv =
_pressure/sqrt_a*(u_gradphi(0)*A_2(2) - A_1(2)*u_gradphi(1));
155 const libMesh::Real dt0_dw =
_pressure/sqrt_a*(A_1(1)*u_gradphi(1) - u_gradphi(0)*A_2(1));
157 const libMesh::Real dt1_du =
_pressure/sqrt_a*(A_1(2)*u_gradphi(1) - u_gradphi(0)*A_2(2));
158 const libMesh::Real dt1_dw =
_pressure/sqrt_a*(u_gradphi(0)*A_2(0) - A_1(0)*u_gradphi(1));
160 const libMesh::Real dt2_du =
_pressure/sqrt_a*(u_gradphi(0)*A_2(1) - A_1(1)*u_gradphi(1));
161 const libMesh::Real dt2_dv =
_pressure/sqrt_a*(A_1(0)*u_gradphi(1) - u_gradphi(0)*A_2(0));
163 Kuv(i,j) -= dt0_dv*u_phi[i][qp]*jac;
164 Kuw(i,j) -= dt0_dw*u_phi[i][qp]*jac;
166 Kvu(i,j) -= dt1_du*u_phi[i][qp]*jac;
167 Kvw(i,j) -= dt1_dw*u_phi[i][qp]*jac;
169 Kwu(i,j) -= dt2_du*u_phi[i][qp]*jac;
170 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.
DisplacementFEVariables _disp_vars
virtual void element_time_derivative(bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
Time dependent part(s) of physics for element interiors.
static PhysicsName elastic_membrane_constant_pressure()
static PhysicsName elastic_membrane()
void reset_pressure(libMesh::Real pressure_in)
void parse_enabled_subdomains(const GetPot &input, const std::string &physics_name)
const libMesh::FEGenericBase< libMesh::Real > * get_fe(const AssemblyContext &context)
ElasticMembraneConstantPressure()