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/fe_interface.h"
45 template<
typename StressStrainLaw>
47 bool is_compressible )
49 _stress_strain_law(input),
51 _is_compressible(is_compressible)
54 if( !input.have_variable(
"Physics/"+physics_name+
"/A") )
56 std::cerr <<
"Error: Must specify initial area for "+physics_name << std::endl
57 <<
" Input the option Physics/"+physics_name+
"/A" << std::endl;
62 (
_A, input,
"Physics/"+physics_name+
"/A",
_A );
71 template<
typename StressStrainLaw>
78 template<
typename StressStrainLaw>
82 std::string section =
"Physics/"+
elastic_cable+
"/output_vars";
84 if( input.have_variable(section) )
86 unsigned int n_vars = input.vector_variable_size(section);
88 for(
unsigned int v = 0; v < n_vars; v++ )
90 std::string name = input(section,
"DIE!",v);
92 if( name == std::string(
"stress") )
96 _stress_indices.resize(1);
101 else if( name == std::string(
"strain") )
105 _strain_indices.resize(1);
110 else if( name == std::string(
"force") )
114 _force_indices.resize(1);
121 std::cerr <<
"Error: Invalue output_vars value for "+
elastic_cable << std::endl
122 <<
" Found " << name << std::endl
123 <<
" Acceptable values are: stress" << std::endl
124 <<
" strain" << std::endl
125 <<
" force " << std::endl;
133 template<
typename StressStrainLaw>
138 const unsigned int n_u_dofs = context.get_dof_indices(_disp_vars.u_var()).size();
140 const std::vector<libMesh::Real> &JxW =
141 this->get_fe(context)->get_JxW();
144 libMesh::DenseSubVector<libMesh::Number> &Fu = context.get_elem_residual(_disp_vars.u_var());
145 libMesh::DenseSubVector<libMesh::Number> &Fv = context.get_elem_residual(_disp_vars.v_var());
146 libMesh::DenseSubVector<libMesh::Number> &Fw = context.get_elem_residual(_disp_vars.w_var());
150 libMesh::DenseSubMatrix<libMesh::Number> &Kuu = context.get_elem_jacobian(_disp_vars.u_var(),_disp_vars.u_var());
151 libMesh::DenseSubMatrix<libMesh::Number> &Kuv = context.get_elem_jacobian(_disp_vars.u_var(),_disp_vars.v_var());
152 libMesh::DenseSubMatrix<libMesh::Number> &Kuw = context.get_elem_jacobian(_disp_vars.u_var(),_disp_vars.w_var());
153 libMesh::DenseSubMatrix<libMesh::Number> &Kvu = context.get_elem_jacobian(_disp_vars.v_var(),_disp_vars.u_var());
154 libMesh::DenseSubMatrix<libMesh::Number> &Kvv = context.get_elem_jacobian(_disp_vars.v_var(),_disp_vars.v_var());
155 libMesh::DenseSubMatrix<libMesh::Number> &Kvw = context.get_elem_jacobian(_disp_vars.v_var(),_disp_vars.w_var());
156 libMesh::DenseSubMatrix<libMesh::Number> &Kwu = context.get_elem_jacobian(_disp_vars.w_var(),_disp_vars.u_var());
157 libMesh::DenseSubMatrix<libMesh::Number> &Kwv = context.get_elem_jacobian(_disp_vars.w_var(),_disp_vars.v_var());
158 libMesh::DenseSubMatrix<libMesh::Number> &Kww = context.get_elem_jacobian(_disp_vars.w_var(),_disp_vars.w_var());
161 unsigned int n_qpoints = context.get_element_qrule().n_points();
164 const std::vector<std::vector<libMesh::Real> >& dphi_dxi = this->get_fe(context)->get_dphidxi();
166 const libMesh::DenseSubVector<libMesh::Number>& u_coeffs = context.get_elem_solution( _disp_vars.u_var() );
167 const libMesh::DenseSubVector<libMesh::Number>& v_coeffs = context.get_elem_solution( _disp_vars.v_var() );
168 const libMesh::DenseSubVector<libMesh::Number>& w_coeffs = context.get_elem_solution( _disp_vars.w_var() );
171 const std::vector<libMesh::RealGradient>& dxdxi = this->get_fe(context)->get_dxyzdxi();
173 const unsigned int dim = 1;
176 for (
unsigned int qp=0; qp != n_qpoints; qp++)
179 libMesh::Gradient grad_u, grad_v, grad_w;
181 for(
unsigned int d = 0; d < n_u_dofs; d++ )
183 libMesh::RealGradient u_gradphi( dphi_dxi[d][qp] );
184 grad_u += u_coeffs(d)*u_gradphi;
185 grad_v += v_coeffs(d)*u_gradphi;
186 grad_w += w_coeffs(d)*u_gradphi;
189 libMesh::RealGradient grad_x( dxdxi[qp](0) );
190 libMesh::RealGradient grad_y( dxdxi[qp](1) );
191 libMesh::RealGradient grad_z( dxdxi[qp](2) );
194 libMesh::Real lambda_sq = 0;
196 this->compute_metric_tensors( qp, *(this->get_fe(context)), context,
197 grad_u, grad_v, grad_w,
198 a_cov, a_contra, A_cov, A_contra,
205 _stress_strain_law.compute_stress_and_elasticity(dim,a_contra,a_cov,A_contra,A_cov,tau,C);
208 libMesh::Real jac = JxW[qp];
210 for (
unsigned int i=0; i != n_u_dofs; i++)
212 libMesh::RealGradient u_gradphi( dphi_dxi[i][qp] );
214 Fu(i) -= tau(0,0)*_A*( (grad_x(0) + grad_u(0))*u_gradphi(0) ) * jac;
216 Fv(i) -= tau(0,0)*_A*( (grad_y(0) + grad_v(0))*u_gradphi(0) ) * jac;
218 Fw(i) -= tau(0,0)*_A*( (grad_z(0) + grad_w(0))*u_gradphi(0) ) * jac;
222 if( compute_jacobian )
224 for(
unsigned int i=0; i != n_u_dofs; i++)
226 libMesh::RealGradient u_gradphi_I( dphi_dxi[i][qp] );
227 for(
unsigned int j=0; j != n_u_dofs; j++)
229 libMesh::RealGradient u_gradphi_J( dphi_dxi[j][qp] );
231 const libMesh::Real diag_term = _A*jac*tau(0,0)*( u_gradphi_J(0)*u_gradphi_I(0));
233 Kuu(i,j) -= diag_term;
235 Kvv(i,j) -= diag_term;
237 Kww(i,j) -= diag_term;
239 const libMesh::Real dgamma_du = ( u_gradphi_J(0)*(grad_x(0)+grad_u(0)) );
241 const libMesh::Real dgamma_dv = ( u_gradphi_J(0)*(grad_y(0)+grad_v(0)) );
243 const libMesh::Real dgamma_dw = ( u_gradphi_J(0)*(grad_z(0)+grad_w(0)) );
245 const libMesh::Real C1 = _A*jac*C(0,0,0,0);
247 const libMesh::Real x_term = C1*( (grad_x(0)+grad_u(0))*u_gradphi_I(0) );
249 const libMesh::Real y_term = C1*( (grad_y(0)+grad_v(0))*u_gradphi_I(0) );
251 const libMesh::Real z_term = C1*( (grad_z(0)+grad_w(0))*u_gradphi_I(0) );
253 Kuu(i,j) -= x_term*dgamma_du;
255 Kuv(i,j) -= x_term*dgamma_dv;
257 Kuw(i,j) -= x_term*dgamma_dw;
259 Kvu(i,j) -= y_term*dgamma_du;
261 Kvv(i,j) -= y_term*dgamma_dv;
263 Kvw(i,j) -= y_term*dgamma_dw;
265 Kwu(i,j) -= z_term*dgamma_du;
267 Kwv(i,j) -= z_term*dgamma_dv;
269 Kww(i,j) -= z_term*dgamma_dw;
280 template<
typename StressStrainLaw>
285 std::vector<BoundaryID> ids = context.side_boundary_ids();
287 for( std::vector<BoundaryID>::const_iterator it = ids.begin();
288 it != ids.end(); it++ )
290 libmesh_assert (*it != libMesh::BoundaryInfo::invalid_id);
292 _bc_handler->apply_neumann_bcs( context, cache, compute_jacobian, *it );
298 template<
typename StressStrainLaw>
303 libmesh_not_implemented();
307 template<
typename StressStrainLaw>
310 const libMesh::Point& point,
311 libMesh::Real& value )
313 bool is_strain =
false;
314 if( !_strain_indices.empty() )
315 is_strain = ( _strain_indices[0] == quantity_index );
317 bool is_stress =
false;
318 if( !_stress_indices.empty() )
319 is_stress = ( _stress_indices[0] == quantity_index );
321 bool is_force =
false;
322 if( !_force_indices.empty() )
323 is_force = ( _force_indices[0] == quantity_index );
325 if( is_strain || is_stress || is_force )
327 const unsigned int n_u_dofs = context.get_dof_indices(_disp_vars.u_var()).size();
329 const libMesh::DenseSubVector<libMesh::Number>& u_coeffs = context.get_elem_solution( _disp_vars.u_var() );
330 const libMesh::DenseSubVector<libMesh::Number>& v_coeffs = context.get_elem_solution( _disp_vars.v_var() );
331 const libMesh::DenseSubVector<libMesh::Number>& w_coeffs = context.get_elem_solution( _disp_vars.w_var() );
334 libMesh::AutoPtr<libMesh::FEGenericBase<libMesh::Real> > fe_new = this->build_new_fe( context.get_elem(), this->get_fe(context), point );
336 const std::vector<std::vector<libMesh::Real> >& dphi_dxi = fe_new->get_dphidxi();
339 const std::vector<libMesh::RealGradient>& dxdxi = fe_new->get_dxyzdxi();
342 libMesh::Gradient grad_u, grad_v, grad_w;
343 for(
unsigned int d = 0; d < n_u_dofs; d++ )
345 libMesh::RealGradient u_gradphi( dphi_dxi[d][0] );
346 grad_u += u_coeffs(d)*u_gradphi;
347 grad_v += v_coeffs(d)*u_gradphi;
348 grad_w += w_coeffs(d)*u_gradphi;
351 libMesh::RealGradient grad_x( dxdxi[0](0) );
352 libMesh::RealGradient grad_y( dxdxi[0](1) );
353 libMesh::RealGradient grad_z( dxdxi[0](2) );
356 libMesh::Real lambda_sq = 0;
358 this->compute_metric_tensors(0, *fe_new, context, grad_u, grad_v, grad_w, a_cov, a_contra, A_cov, A_contra, lambda_sq );
360 libMesh::Real det_a = a_cov(0,0)*a_cov(1,1) - a_cov(0,1)*a_cov(1,0);
361 libMesh::Real det_A = A_cov(0,0)*A_cov(1,1) - A_cov(0,1)*A_cov(1,0);
363 libMesh::Real I3 = lambda_sq*det_A/det_a;
366 _stress_strain_law.compute_stress(1,a_contra,a_cov,A_contra,A_cov,tau);
371 if( _strain_indices[0] == quantity_index )
373 value = 0.5*(A_cov(0,0) - a_cov(0,0));
385 if( _stress_indices[0] == quantity_index )
388 value = tau(0,0)/std::sqrt(I3);
398 if( _force_indices[0] == quantity_index )
401 value = tau(0,0)/std::sqrt(I3)*_A;
413 template<
typename StressStrainLaw>
415 const libMesh::FEGenericBase<libMesh::Real>* fe,
416 const libMesh::Point p )
419 FEType fe_type = fe->get_fe_type();
424 libmesh_assert(&elem || fe_type.family == SCALAR);
426 unsigned int elem_dim = &elem ? elem.dim() : 0;
428 AutoPtr<FEGenericBase<libMesh::Real> >
429 fe_new(FEGenericBase<libMesh::Real>::build(elem_dim, fe_type));
433 Point master_point = &elem ?
434 FEInterface::inverse_map(elem_dim, fe_type, &elem, p) :
437 std::vector<Point> coor(1, master_point);
440 fe_new->reinit (&elem, &coor);
445 template<
typename StressStrainLaw>
447 const libMesh::FEBase& elem,
449 const libMesh::Gradient& grad_u,
450 const libMesh::Gradient& grad_v,
451 const libMesh::Gradient& grad_w,
456 libMesh::Real& lambda_sq )
458 const std::vector<libMesh::RealGradient>& dxdxi = elem.get_dxyzdxi();
460 const std::vector<libMesh::Real>& dxidx = elem.get_dxidx();
461 const std::vector<libMesh::Real>& dxidy = elem.get_dxidy();
462 const std::vector<libMesh::Real>& dxidz = elem.get_dxidz();
464 libMesh::RealGradient dxi( dxidx[qp], dxidy[qp], dxidz[qp] );
466 libMesh::RealGradient dudxi( grad_u(0), grad_v(0), grad_w(0) );
470 a_cov(0,0) = dxdxi[qp]*dxdxi[qp];
476 A_cov(0,0) = (dxdxi[qp] + dudxi)*(dxdxi[qp] + dudxi);
480 a_contra(0,0) = 1/a_cov(0,0);
486 A_contra(0,0) = 1/A_cov(0,0);
489 if( _is_compressible )
491 libmesh_not_implemented();
497 lambda_sq = a_cov(0,0)/A_cov(0,0);
501 A_cov(1,1) = lambda_sq;
502 A_cov(2,2) = lambda_sq;
503 A_contra(1,1) = 1.0/lambda_sq;
504 A_contra(2,2) = 1.0/lambda_sq;
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.
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)
unsigned int register_quantity(std::string name)
Register quantity to be postprocessed.
GRINS::ICHandlingBase * _ic_handler
GRINS::BCHandlingBase * _bc_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.
Base class for reading and handling initial conditions for physics classes.
virtual void side_time_derivative(bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
Time dependent part(s) of physics for boundaries of elements on the domain boundary.
virtual void register_postprocessing_vars(const GetPot &input, PostProcessedQuantities< libMesh::Real > &postprocessing)
Register postprocessing variables for ElasticCable.
const PhysicsName elastic_cable
virtual void mass_residual(bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
Mass matrix part(s) for element interiors. All boundary terms lie within the time_derivative part...
virtual void element_time_derivative(bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
Time dependent part(s) of physics for element interiors.
libMesh::AutoPtr< libMesh::FEGenericBase< libMesh::Real > > build_new_fe(const libMesh::Elem &elem, const libMesh::FEGenericBase< libMesh::Real > *fe, const libMesh::Point p)