29 #include "grins_config.h"
37 #include "libmesh/getpot.h"
38 #include "libmesh/quadrature.h"
39 #include "libmesh/boundary_info.h"
40 #include "libmesh/fem_system.h"
41 #include "libmesh/elem.h"
45 template<
typename StressStrainLaw>
47 bool is_compressible )
53 template<
typename StressStrainLaw>
59 if( input.have_variable(section) )
61 unsigned int n_vars = input.vector_variable_size(section);
63 for(
unsigned int v = 0; v < n_vars; v++ )
65 std::string name = input(section,
"DIE!",v);
67 if( name == std::string(
"stress") )
71 _stress_indices.resize(7);
87 else if( name == std::string(
"stress_zz" ) )
92 else if( name == std::string(
"strain") )
95 _strain_indices.resize(3);
106 <<
" Found " << name << std::endl
107 <<
" Acceptable values are: stress" << std::endl
108 <<
" strain" << std::endl;
115 template<
typename StressStrainLaw>
120 const unsigned int n_u_dofs = context.get_dof_indices(this->_disp_vars.u()).size();
122 const std::vector<libMesh::Real> &JxW =
123 this->get_fe(context)->get_JxW();
126 libMesh::DenseSubVector<libMesh::Number> &Fu = context.get_elem_residual(this->_disp_vars.u());
127 libMesh::DenseSubVector<libMesh::Number> &Fv = context.get_elem_residual(this->_disp_vars.v());
128 libMesh::DenseSubVector<libMesh::Number> &Fw = context.get_elem_residual(this->_disp_vars.w());
130 libMesh::DenseSubMatrix<libMesh::Number>& Kuu = context.get_elem_jacobian(this->_disp_vars.u(),this->_disp_vars.u());
131 libMesh::DenseSubMatrix<libMesh::Number>& Kuv = context.get_elem_jacobian(this->_disp_vars.u(),this->_disp_vars.v());
132 libMesh::DenseSubMatrix<libMesh::Number>& Kuw = context.get_elem_jacobian(this->_disp_vars.u(),this->_disp_vars.w());
134 libMesh::DenseSubMatrix<libMesh::Number>& Kvu = context.get_elem_jacobian(this->_disp_vars.v(),this->_disp_vars.u());
135 libMesh::DenseSubMatrix<libMesh::Number>& Kvv = context.get_elem_jacobian(this->_disp_vars.v(),this->_disp_vars.v());
136 libMesh::DenseSubMatrix<libMesh::Number>& Kvw = context.get_elem_jacobian(this->_disp_vars.v(),this->_disp_vars.w());
138 libMesh::DenseSubMatrix<libMesh::Number>& Kwu = context.get_elem_jacobian(this->_disp_vars.w(),this->_disp_vars.u());
139 libMesh::DenseSubMatrix<libMesh::Number>& Kwv = context.get_elem_jacobian(this->_disp_vars.w(),this->_disp_vars.v());
140 libMesh::DenseSubMatrix<libMesh::Number>& Kww = context.get_elem_jacobian(this->_disp_vars.w(),this->_disp_vars.w());
142 unsigned int n_qpoints = context.get_element_qrule().n_points();
145 const std::vector<std::vector<libMesh::Real> >& dphi_dxi =
146 this->get_fe(context)->get_dphidxi();
148 const std::vector<std::vector<libMesh::Real> >& dphi_deta =
149 this->get_fe(context)->get_dphideta();
151 const libMesh::DenseSubVector<libMesh::Number>& u_coeffs = context.get_elem_solution( this->_disp_vars.u() );
152 const libMesh::DenseSubVector<libMesh::Number>& v_coeffs = context.get_elem_solution( this->_disp_vars.v() );
153 const libMesh::DenseSubVector<libMesh::Number>& w_coeffs = context.get_elem_solution( this->_disp_vars.w() );
156 const std::vector<libMesh::RealGradient>& dxdxi = this->get_fe(context)->get_dxyzdxi();
157 const std::vector<libMesh::RealGradient>& dxdeta = this->get_fe(context)->get_dxyzdeta();
159 const unsigned int dim = 2;
161 for (
unsigned int qp=0; qp != n_qpoints; qp++)
164 libMesh::Gradient grad_u, grad_v, grad_w;
165 for(
unsigned int d = 0; d < n_u_dofs; d++ )
167 libMesh::RealGradient u_gradphi( dphi_dxi[d][qp], dphi_deta[d][qp] );
168 grad_u += u_coeffs(d)*u_gradphi;
169 grad_v += v_coeffs(d)*u_gradphi;
170 grad_w += w_coeffs(d)*u_gradphi;
173 libMesh::RealGradient grad_x( dxdxi[qp](0), dxdeta[qp](0) );
174 libMesh::RealGradient grad_y( dxdxi[qp](1), dxdeta[qp](1) );
175 libMesh::RealGradient grad_z( dxdxi[qp](2), dxdeta[qp](2) );
179 libMesh::Real lambda_sq = 0;
181 this->compute_metric_tensors( qp, *(this->get_fe(context)), context,
182 grad_u, grad_v, grad_w,
183 a_cov, a_contra, A_cov, A_contra,
189 this->_stress_strain_law.compute_stress_and_elasticity(dim,a_contra,a_cov,A_contra,A_cov,tau,C);
191 libMesh::Real jac = JxW[qp];
193 for (
unsigned int i=0; i != n_u_dofs; i++)
195 libMesh::RealGradient u_gradphi( dphi_dxi[i][qp], dphi_deta[i][qp] );
197 for(
unsigned int alpha = 0; alpha < dim; alpha++ )
199 for(
unsigned int beta = 0; beta < dim; beta++ )
201 libMesh::Real factor = 0.5*tau(alpha,beta)*this->_h0*jac;
203 Fu(i) += factor*( (grad_x(beta) + grad_u(beta))*u_gradphi(alpha) +
204 (grad_x(alpha) + grad_u(alpha))*u_gradphi(beta) );
206 Fv(i) += factor*( (grad_y(beta) + grad_v(beta))*u_gradphi(alpha) +
207 (grad_y(alpha) + grad_v(alpha))*u_gradphi(beta) );
209 Fw(i) += factor*( (grad_z(beta) + grad_w(beta))*u_gradphi(alpha) +
210 (grad_z(alpha) + grad_w(alpha))*u_gradphi(beta) );
215 if( compute_jacobian )
217 for (
unsigned int i=0; i != n_u_dofs; i++)
219 libMesh::RealGradient u_gradphi_i( dphi_dxi[i][qp], dphi_deta[i][qp] );
221 for (
unsigned int j=0; j != n_u_dofs; j++)
223 libMesh::RealGradient u_gradphi_j( dphi_dxi[j][qp], dphi_deta[j][qp] );
225 for(
unsigned int alpha = 0; alpha < dim; alpha++ )
227 for(
unsigned int beta = 0; beta < dim; beta++ )
229 const libMesh::Real diag_term = 0.5*this->_h0*jac*tau(alpha,beta)*context.get_elem_solution_derivative()*
230 ( u_gradphi_j(beta)*u_gradphi_i(alpha) +
231 u_gradphi_j(alpha)*u_gradphi_i(beta) );
232 Kuu(i,j) += diag_term;
234 Kvv(i,j) += diag_term;
236 Kww(i,j) += diag_term;
238 for(
unsigned int lambda = 0; lambda < dim; lambda++ )
240 for(
unsigned int mu = 0; mu < dim; mu++ )
242 const libMesh::Real dgamma_du = 0.5*( u_gradphi_j(lambda)*(grad_x(mu)+grad_u(mu)) +
243 (grad_x(lambda)+grad_u(lambda))*u_gradphi_j(mu) );
245 const libMesh::Real dgamma_dv = 0.5*( u_gradphi_j(lambda)*(grad_y(mu)+grad_v(mu)) +
246 (grad_y(lambda)+grad_v(lambda))*u_gradphi_j(mu) );
248 const libMesh::Real dgamma_dw = 0.5*( u_gradphi_j(lambda)*(grad_z(mu)+grad_w(mu)) +
249 (grad_z(lambda)+grad_w(lambda))*u_gradphi_j(mu) );
251 const libMesh::Real C1 = 0.5*this->_h0*jac*C(alpha,beta,lambda,mu)*context.get_elem_solution_derivative();
253 const libMesh::Real x_term = C1*( (grad_x(beta)+grad_u(beta))*u_gradphi_i(alpha) +
254 (grad_x(alpha)+grad_u(alpha))*u_gradphi_i(beta) );
256 const libMesh::Real y_term = C1*( (grad_y(beta)+grad_v(beta))*u_gradphi_i(alpha) +
257 (grad_y(alpha)+grad_v(alpha))*u_gradphi_i(beta) );
259 const libMesh::Real z_term = C1*( (grad_z(beta)+grad_w(beta))*u_gradphi_i(alpha) +
260 (grad_z(alpha)+grad_w(alpha))*u_gradphi_i(beta) );
262 Kuu(i,j) += x_term*dgamma_du;
264 Kuv(i,j) += x_term*dgamma_dv;
266 Kuw(i,j) += x_term*dgamma_dw;
268 Kvu(i,j) += y_term*dgamma_du;
270 Kvv(i,j) += y_term*dgamma_dv;
272 Kvw(i,j) += y_term*dgamma_dw;
274 Kwu(i,j) += z_term*dgamma_du;
276 Kwv(i,j) += z_term*dgamma_dv;
278 Kww(i,j) += z_term*dgamma_dw;
290 template<
typename StressStrainLaw>
296 if( this->_is_compressible )
298 unsigned int n_qpoints = context.get_element_qrule().n_points();
300 const unsigned int n_u_dofs = context.get_dof_indices(this->_disp_vars.u()).size();
302 const std::vector<libMesh::Real> &JxW = context.get_element_fe(this->_lambda_sq_var)->get_JxW();
304 libMesh::DenseSubVector<libMesh::Number>& Fl = context.get_elem_residual(this->_lambda_sq_var);
306 const std::vector<std::vector<libMesh::Real> >& phi =
307 context.get_element_fe(this->_lambda_sq_var)->get_phi();
309 const unsigned int n_lambda_sq_dofs = context.get_dof_indices(this->_lambda_sq_var).size();
311 const libMesh::DenseSubVector<libMesh::Number>& u_coeffs = context.get_elem_solution( this->_disp_vars.u() );
312 const libMesh::DenseSubVector<libMesh::Number>& v_coeffs = context.get_elem_solution( this->_disp_vars.v() );
313 const libMesh::DenseSubVector<libMesh::Number>& w_coeffs = context.get_elem_solution( this->_disp_vars.w() );
316 const std::vector<std::vector<libMesh::Real> >& dphi_dxi =
317 this->get_fe(context)->get_dphidxi();
319 const std::vector<std::vector<libMesh::Real> >& dphi_deta =
320 this->get_fe(context)->get_dphideta();
322 for (
unsigned int qp=0; qp != n_qpoints; qp++)
324 libMesh::Real jac = JxW[qp];
326 libMesh::Gradient grad_u, grad_v, grad_w;
327 for(
unsigned int d = 0; d < n_u_dofs; d++ )
329 libMesh::RealGradient u_gradphi( dphi_dxi[d][qp], dphi_deta[d][qp] );
330 grad_u += u_coeffs(d)*u_gradphi;
331 grad_v += v_coeffs(d)*u_gradphi;
332 grad_w += w_coeffs(d)*u_gradphi;
336 libMesh::Real lambda_sq = 0;
338 this->compute_metric_tensors( qp, *(this->get_fe(context)), context,
339 grad_u, grad_v, grad_w,
340 a_cov, a_contra, A_cov, A_contra,
343 libMesh::Real stress_33 = this->_stress_strain_law.compute_33_stress( a_contra, a_cov, A_contra, A_cov );
345 for (
unsigned int i=0; i != n_lambda_sq_dofs; i++)
347 Fl(i) += stress_33*phi[i][qp]*jac;
350 if( compute_jacobian )
352 libmesh_not_implemented();
358 template<
typename StressStrainLaw>
361 const libMesh::Point& point,
362 libMesh::Real& value )
364 bool is_stress =
false;
365 if( !_stress_indices.empty() )
366 is_stress= ( _stress_indices[0] == quantity_index ||
367 _stress_indices[1] == quantity_index ||
368 _stress_indices[2] == quantity_index ||
369 _stress_indices[3] == quantity_index ||
370 _stress_indices[4] == quantity_index ||
371 _stress_indices[5] == quantity_index ||
372 _stress_indices[6] == quantity_index ||
373 _stress_zz_index == quantity_index );
375 bool is_strain =
false;
376 if( !_strain_indices.empty() )
377 is_strain = ( _strain_indices[0] == quantity_index ||
378 _strain_indices[1] == quantity_index ||
379 _strain_indices[2] == quantity_index );
381 if( is_stress || is_strain )
383 const unsigned int n_u_dofs = context.get_dof_indices(this->_disp_vars.u()).size();
385 const libMesh::DenseSubVector<libMesh::Number>& u_coeffs = context.get_elem_solution( this->_disp_vars.u() );
386 const libMesh::DenseSubVector<libMesh::Number>& v_coeffs = context.get_elem_solution( this->_disp_vars.v() );
387 const libMesh::DenseSubVector<libMesh::Number>& w_coeffs = context.get_elem_solution( this->_disp_vars.w() );
390 libMesh::UniquePtr<libMesh::FEGenericBase<libMesh::Real> > fe_new =
391 this->build_new_fe( &context.get_elem(), this->get_fe(context),
394 const std::vector<std::vector<libMesh::Real> >& dphi_dxi =
395 fe_new->get_dphidxi();
397 const std::vector<std::vector<libMesh::Real> >& dphi_deta =
398 fe_new->get_dphideta();
401 const std::vector<libMesh::RealGradient>& dxdxi = fe_new->get_dxyzdxi();
402 const std::vector<libMesh::RealGradient>& dxdeta = fe_new->get_dxyzdeta();
405 libMesh::Gradient grad_u, grad_v, grad_w;
406 for(
unsigned int d = 0; d < n_u_dofs; d++ )
408 libMesh::RealGradient u_gradphi( dphi_dxi[d][0], dphi_deta[d][0] );
409 grad_u += u_coeffs(d)*u_gradphi;
410 grad_v += v_coeffs(d)*u_gradphi;
411 grad_w += w_coeffs(d)*u_gradphi;
414 libMesh::RealGradient grad_x( dxdxi[0](0), dxdeta[0](0) );
415 libMesh::RealGradient grad_y( dxdxi[0](1), dxdeta[0](1) );
416 libMesh::RealGradient grad_z( dxdxi[0](2), dxdeta[0](2) );
419 libMesh::Real lambda_sq = 0;
422 this->compute_metric_tensors(0, *fe_new, context, grad_u, grad_v, grad_w,
423 a_cov, a_contra, A_cov, A_contra, lambda_sq );
428 if( _strain_indices[0] == quantity_index )
430 value = 0.5*(A_cov(0,0) - a_cov(0,0));
432 else if( _strain_indices[1] == quantity_index )
434 value = 0.5*(A_cov(0,1) - a_cov(0,1));
436 else if( _strain_indices[2] == quantity_index )
438 value = 0.5*(A_cov(1,1) - a_cov(1,1));
450 libMesh::Real det_a = a_cov(0,0)*a_cov(1,1) - a_cov(0,1)*a_cov(1,0);
451 libMesh::Real det_A = A_cov(0,0)*A_cov(1,1) - A_cov(0,1)*A_cov(1,0);
453 libMesh::Real I3 = lambda_sq*det_A/det_a;
456 this->_stress_strain_law.compute_stress(2,a_contra,a_cov,A_contra,A_cov,tau);
458 if( _stress_indices[0] == quantity_index )
461 value = tau(0,0)/std::sqrt(I3);
463 else if( _stress_indices[1] == quantity_index )
466 value = tau(0,1)/std::sqrt(I3);
468 else if( _stress_indices[2] == quantity_index )
471 value = tau(1,1)/std::sqrt(I3);
473 else if( _stress_indices[3] == quantity_index )
475 value = 0.5*(tau(0,0) + tau(1,1)) + std::sqrt(0.25*(tau(0,0)-tau(1,1))*(tau(0,0)-tau(1,1))
476 + tau(0,1)*tau(0,1) );
478 else if( _stress_indices[4] == quantity_index ||
479 _stress_indices[5] == quantity_index ||
480 _stress_indices[6] == quantity_index )
482 if(this->_is_compressible)
484 tau(2,2) = this->_stress_strain_law.compute_33_stress( a_contra, a_cov, A_contra, A_cov );
487 libMesh::Real stress_I1 = tau(0,0) + tau(1,1) + tau(2,2);
488 libMesh::Real stress_I2 = 0.5*(stress_I1*stress_I1 - (tau(0,0)*tau(0,0) + tau(1,1)*tau(1,1)
489 + tau(2,2)*tau(2,2) + tau(0,1)*tau(0,1)
490 + tau(1,0)*tau(1,0)) );
492 libMesh::Real stress_I3 = tau(2,2)*( tau(0,0)*tau(1,1) - tau(1,0)*tau(0,1) );
498 libMesh::Real C1 = (stress_I1*stress_I1 - 3*stress_I2);
501 libMesh::Real C2 = (2*stress_I1*stress_I1*stress_I1 - 9*stress_I1*stress_I2 + 27*stress_I3)/54;
503 libMesh::Real theta = std::acos( C2/(2*std::sqrt(C1*C1*C1)) )/3.0;
505 if( _stress_indices[4] == quantity_index )
507 value = (stress_I1 + 2.0*std::sqrt(C1)*std::cos(theta))/3.0;
510 if( _stress_indices[5] == quantity_index )
512 value = (stress_I1 + 2.0*std::sqrt(C1)*std::cos(theta+
Constants::two_pi/3.0))/3.0;
515 if( _stress_indices[6] == quantity_index )
517 value = (stress_I1 + 2.0*std::sqrt(C1)*std::cos(theta+2.0*
Constants::two_pi/3.0))/3.0;
520 else if( _stress_zz_index == quantity_index )
522 value = this->_stress_strain_law.compute_33_stress( a_contra, a_cov, A_contra, A_cov );
unsigned int register_quantity(std::string name)
Register quantity to be postprocessed.
GRINS::ICHandlingBase * _ic_handler
virtual void element_time_derivative(bool compute_jacobian, AssemblyContext &context, CachedValues &)
Time dependent part(s) of physics for element interiors.
Base class for reading and handling initial conditions for physics classes.
virtual void element_constraint(bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
Constraint part(s) of physics for element interiors.
static PhysicsName elastic_membrane()
virtual void register_postprocessing_vars(const GetPot &input, PostProcessedQuantities< libMesh::Real > &postprocessing)
Register postprocessing variables for ElasticMembrane.
virtual void compute_postprocessed_quantity(unsigned int quantity_index, const AssemblyContext &context, const libMesh::Point &point, libMesh::Real &value)
Compute the registered postprocessed quantities.
const libMesh::Real two_pi