29 #include "grins_config.h"
35 #include "libmesh/getpot.h"
36 #include "libmesh/quadrature.h"
37 #include "libmesh/boundary_info.h"
38 #include "libmesh/fem_system.h"
39 #include "libmesh/elem.h"
43 template<
typename StressStrainLaw>
45 bool is_compressible )
51 template<
typename StressStrainLaw>
57 if( input.have_variable(section) )
59 unsigned int n_vars = input.vector_variable_size(section);
61 for(
unsigned int v = 0; v < n_vars; v++ )
63 std::string name = input(section,
"DIE!",v);
65 if( name == std::string(
"stress") )
69 _stress_indices.resize(1);
74 else if( name == std::string(
"strain") )
78 _strain_indices.resize(1);
83 else if( name == std::string(
"force") )
87 _force_indices.resize(1);
95 <<
" Found " << name << std::endl
96 <<
" Acceptable values are: stress" << std::endl
97 <<
" strain" << std::endl
98 <<
" force " << std::endl;
106 template<
typename StressStrainLaw>
109 const libMesh::Point& point,
110 libMesh::Real& value )
112 bool is_strain =
false;
113 if( !_strain_indices.empty() )
114 is_strain = ( _strain_indices[0] == quantity_index );
116 bool is_stress =
false;
117 if( !_stress_indices.empty() )
118 is_stress = ( _stress_indices[0] == quantity_index );
120 bool is_force =
false;
121 if( !_force_indices.empty() )
122 is_force = ( _force_indices[0] == quantity_index );
124 if( is_strain || is_stress || is_force )
126 const unsigned int n_u_dofs = context.get_dof_indices(this->_disp_vars.u()).size();
128 const libMesh::DenseSubVector<libMesh::Number>& u_coeffs = context.get_elem_solution( this->_disp_vars.u() );
129 const libMesh::DenseSubVector<libMesh::Number>& v_coeffs = context.get_elem_solution( this->_disp_vars.v() );
130 const libMesh::DenseSubVector<libMesh::Number>& w_coeffs = context.get_elem_solution( this->_disp_vars.w() );
133 libMesh::UniquePtr<libMesh::FEGenericBase<libMesh::Real> > fe_new = this->build_new_fe( &context.get_elem(), this->get_fe(context), point );
135 const std::vector<std::vector<libMesh::Real> >& dphi_dxi = fe_new->get_dphidxi();
138 const std::vector<libMesh::RealGradient>& dxdxi = fe_new->get_dxyzdxi();
141 libMesh::Gradient grad_u, grad_v, grad_w;
142 for(
unsigned int d = 0; d < n_u_dofs; d++ )
144 libMesh::RealGradient u_gradphi( dphi_dxi[d][0] );
145 grad_u += u_coeffs(d)*u_gradphi;
146 grad_v += v_coeffs(d)*u_gradphi;
147 grad_w += w_coeffs(d)*u_gradphi;
150 libMesh::RealGradient grad_x( dxdxi[0](0) );
151 libMesh::RealGradient grad_y( dxdxi[0](1) );
152 libMesh::RealGradient grad_z( dxdxi[0](2) );
155 libMesh::Real lambda_sq = 0;
157 this->compute_metric_tensors(0, *fe_new, context, grad_u, grad_v, grad_w, a_cov, a_contra, A_cov, A_contra, lambda_sq );
159 libMesh::Real det_a = a_cov(0,0)*a_cov(1,1) - a_cov(0,1)*a_cov(1,0);
160 libMesh::Real det_A = A_cov(0,0)*A_cov(1,1) - A_cov(0,1)*A_cov(1,0);
162 libMesh::Real I3 = lambda_sq*det_A/det_a;
165 this->_stress_strain_law.compute_stress(1,a_contra,a_cov,A_contra,A_cov,tau);
170 if( _strain_indices[0] == quantity_index )
172 value = 0.5*(A_cov(0,0) - a_cov(0,0));
184 if( _stress_indices[0] == quantity_index )
187 value = tau(0,0)/std::sqrt(I3);
197 if( _force_indices[0] == quantity_index )
200 value = tau(0,0)/std::sqrt(I3)*this->_A;
210 template<
typename StressStrainLaw>
215 const unsigned int n_u_dofs = context.get_dof_indices(this->_disp_vars.u()).size();
217 const std::vector<libMesh::Real> &JxW =
218 this->get_fe(context)->get_JxW();
221 libMesh::DenseSubVector<libMesh::Number> &Fu = context.get_elem_residual(this->_disp_vars.u());
222 libMesh::DenseSubVector<libMesh::Number> &Fv = context.get_elem_residual(this->_disp_vars.v());
223 libMesh::DenseSubVector<libMesh::Number> &Fw = context.get_elem_residual(this->_disp_vars.w());
227 libMesh::DenseSubMatrix<libMesh::Number> &Kuu = context.get_elem_jacobian(this->_disp_vars.u(),this->_disp_vars.u());
228 libMesh::DenseSubMatrix<libMesh::Number> &Kuv = context.get_elem_jacobian(this->_disp_vars.u(),this->_disp_vars.v());
229 libMesh::DenseSubMatrix<libMesh::Number> &Kuw = context.get_elem_jacobian(this->_disp_vars.u(),this->_disp_vars.w());
230 libMesh::DenseSubMatrix<libMesh::Number> &Kvu = context.get_elem_jacobian(this->_disp_vars.v(),this->_disp_vars.u());
231 libMesh::DenseSubMatrix<libMesh::Number> &Kvv = context.get_elem_jacobian(this->_disp_vars.v(),this->_disp_vars.v());
232 libMesh::DenseSubMatrix<libMesh::Number> &Kvw = context.get_elem_jacobian(this->_disp_vars.v(),this->_disp_vars.w());
233 libMesh::DenseSubMatrix<libMesh::Number> &Kwu = context.get_elem_jacobian(this->_disp_vars.w(),this->_disp_vars.u());
234 libMesh::DenseSubMatrix<libMesh::Number> &Kwv = context.get_elem_jacobian(this->_disp_vars.w(),this->_disp_vars.v());
235 libMesh::DenseSubMatrix<libMesh::Number> &Kww = context.get_elem_jacobian(this->_disp_vars.w(),this->_disp_vars.w());
238 unsigned int n_qpoints = context.get_element_qrule().n_points();
241 const std::vector<std::vector<libMesh::Real> >& dphi_dxi = this->get_fe(context)->get_dphidxi();
243 const libMesh::DenseSubVector<libMesh::Number>& u_coeffs = context.get_elem_solution( this->_disp_vars.u() );
244 const libMesh::DenseSubVector<libMesh::Number>& v_coeffs = context.get_elem_solution( this->_disp_vars.v() );
245 const libMesh::DenseSubVector<libMesh::Number>& w_coeffs = context.get_elem_solution( this->_disp_vars.w() );
248 const std::vector<libMesh::RealGradient>& dxdxi = this->get_fe(context)->get_dxyzdxi();
250 const unsigned int dim = 1;
252 for (
unsigned int qp=0; qp != n_qpoints; qp++)
255 libMesh::Gradient grad_u, grad_v, grad_w;
257 for(
unsigned int d = 0; d < n_u_dofs; d++ )
259 libMesh::RealGradient u_gradphi( dphi_dxi[d][qp] );
260 grad_u += u_coeffs(d)*u_gradphi;
261 grad_v += v_coeffs(d)*u_gradphi;
262 grad_w += w_coeffs(d)*u_gradphi;
265 libMesh::RealGradient grad_x( dxdxi[qp](0) );
266 libMesh::RealGradient grad_y( dxdxi[qp](1) );
267 libMesh::RealGradient grad_z( dxdxi[qp](2) );
270 libMesh::Real lambda_sq = 0;
272 this->compute_metric_tensors( qp, *(this->get_fe(context)), context,
273 grad_u, grad_v, grad_w,
274 a_cov, a_contra, A_cov, A_contra,
280 this->_stress_strain_law.compute_stress_and_elasticity(dim,a_contra,a_cov,A_contra,A_cov,tau,C);
283 libMesh::Real jac = JxW[qp];
285 for (
unsigned int i=0; i != n_u_dofs; i++)
287 libMesh::RealGradient u_gradphi( dphi_dxi[i][qp] );
289 const libMesh::Real res_term = tau(0,0)*this->_A*jac*u_gradphi(0);
291 Fu(i) += res_term*(grad_x(0) + grad_u(0));
293 Fv(i) += res_term*(grad_y(0) + grad_v(0));
295 Fw(i) += res_term*(grad_z(0) + grad_w(0));
298 if( compute_jacobian )
300 for(
unsigned int i=0; i != n_u_dofs; i++)
302 libMesh::RealGradient u_gradphi_I( dphi_dxi[i][qp] );
303 for(
unsigned int j=0; j != n_u_dofs; j++)
305 libMesh::RealGradient u_gradphi_J( dphi_dxi[j][qp] );
307 const libMesh::Real diag_term = this->_A*jac*tau(0,0)*( u_gradphi_J(0)*u_gradphi_I(0))*context.get_elem_solution_derivative();
309 Kuu(i,j) += diag_term;
311 Kvv(i,j) += diag_term;
313 Kww(i,j) += diag_term;
315 const libMesh::Real dgamma_du = ( u_gradphi_J(0)*(grad_x(0)+grad_u(0)) );
317 const libMesh::Real dgamma_dv = ( u_gradphi_J(0)*(grad_y(0)+grad_v(0)) );
319 const libMesh::Real dgamma_dw = ( u_gradphi_J(0)*(grad_z(0)+grad_w(0)) );
321 const libMesh::Real C1 = this->_A*jac*C(0,0,0,0)*context.get_elem_solution_derivative();
323 const libMesh::Real x_term = C1*( (grad_x(0)+grad_u(0))*u_gradphi_I(0) );
325 const libMesh::Real y_term = C1*( (grad_y(0)+grad_v(0))*u_gradphi_I(0) );
327 const libMesh::Real z_term = C1*( (grad_z(0)+grad_w(0))*u_gradphi_I(0) );
329 Kuu(i,j) += x_term*dgamma_du;
331 Kuv(i,j) += x_term*dgamma_dv;
333 Kuw(i,j) += x_term*dgamma_dw;
335 Kvu(i,j) += y_term*dgamma_du;
337 Kvv(i,j) += y_term*dgamma_dv;
339 Kvw(i,j) += y_term*dgamma_dw;
341 Kwu(i,j) += z_term*dgamma_du;
343 Kwv(i,j) += z_term*dgamma_dv;
345 Kww(i,j) += z_term*dgamma_dw;
unsigned int register_quantity(std::string name)
Register quantity to be postprocessed.
GRINS::ICHandlingBase * _ic_handler
virtual void compute_postprocessed_quantity(unsigned int quantity_index, const AssemblyContext &context, const libMesh::Point &point, libMesh::Real &value)
Compute the registered postprocessed quantities.
static PhysicsName elastic_cable()
Base class for reading and handling initial conditions for physics classes.
virtual void register_postprocessing_vars(const GetPot &input, PostProcessedQuantities< libMesh::Real > &postprocessing)
Register postprocessing variables for ElasticCable.
virtual void element_time_derivative(bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
Time dependent part(s) of physics for element interiors.