From 32a648c68af97c0f26882b8218701ccd12a16730 Mon Sep 17 00:00:00 2001 From: Jean-Paul Pelteret Date: Tue, 28 Nov 2017 10:47:10 +0200 Subject: [PATCH] Make cook membrane problem dimension independent --- .../cook_membrane.cc | 26 +++++++++---------- 1 file changed, 12 insertions(+), 14 deletions(-) diff --git a/Quasi_static_Finite_strain_Compressible_Elasticity/cook_membrane.cc b/Quasi_static_Finite_strain_Compressible_Elasticity/cook_membrane.cc index 0721148..55a75a6 100644 --- a/Quasi_static_Finite_strain_Compressible_Elasticity/cook_membrane.cc +++ b/Quasi_static_Finite_strain_Compressible_Elasticity/cook_membrane.cc @@ -542,7 +542,7 @@ namespace Cook_Membrane void update_material_data(const Tensor<2, dim> &F) { det_F = determinant(F); - b_bar = std::pow(det_F, -2.0 / 3.0) * symmetrize(F * transpose(F)); + b_bar = std::pow(det_F, -2.0 / dim) * symmetrize(F * transpose(F)); Assert(det_F > 0, ExcInternalError()); } @@ -651,9 +651,9 @@ namespace Cook_Membrane tau_iso); const SymmetricTensor<4, dim> c_bar = get_c_bar(); - return (2.0 / 3.0) * trace(tau_bar) + return (2.0 / dim) * trace(tau_bar) * StandardTensors::dev_P - - (2.0 / 3.0) * (tau_iso_x_I + I_x_tau_iso) + - (2.0 / dim) * (tau_iso_x_I + I_x_tau_iso) + StandardTensors::dev_P * c_bar * StandardTensors::dev_P; } @@ -1735,7 +1735,9 @@ Point grid_y_transform (const Point &pt_in) std::cout << "_"; std::cout << std::endl; - const Point soln_pt (48.0*parameters.scale,60.0*parameters.scale,0.5*parameters.scale); + Point soln_pt (48.0*parameters.scale,60.0*parameters.scale); + if (dim == 3) + soln_pt[2] = 0.5*parameters.scale; double vertical_tip_displacement = 0.0; double vertical_tip_displacement_check = 0.0; @@ -2102,7 +2104,8 @@ Point grid_y_transform (const Point &pt_in) // vector. const double time_ramp = (time.current() / time.end()); const double magnitude = (1.0/(16.0*parameters.scale*1.0*parameters.scale))*time_ramp; // (Total force) / (RHS surface area) - static const Tensor<1, dim> dir ({0.0,1.0,0.0}); + Tensor<1,dim> dir; + dir[1] = 1.0; const Tensor<1, dim> traction = magnitude*dir; for (unsigned int i = 0; i < dofs_per_cell; ++i) @@ -2165,9 +2168,6 @@ Point grid_y_transform (const Point &pt_in) // provide it with an extractor object for the component we wish to // select. To this end we first set up such extractor objects and later // use it when generating the relevant component masks: - const FEValuesExtractors::Scalar x_displacement(0); - const FEValuesExtractors::Scalar y_displacement(1); - const FEValuesExtractors::Scalar z_displacement(2); // Fixed left hand side of the beam { @@ -2178,23 +2178,21 @@ Point grid_y_transform (const Point &pt_in) boundary_id, ZeroFunction(n_components), constraints, - fe.component_mask(x_displacement) | - fe.component_mask(y_displacement) | - fe.component_mask(z_displacement)); + fe.component_mask(u_fe)); else VectorTools::interpolate_boundary_values(dof_handler_ref, boundary_id, ZeroFunction(n_components), constraints, - fe.component_mask(x_displacement) | - fe.component_mask(y_displacement) | - fe.component_mask(z_displacement)); + fe.component_mask(u_fe)); } // Zero Z-displacement through thickness direction // This corresponds to a plane strain condition being imposed on the beam + if (dim == 3) { const int boundary_id = 2; + const FEValuesExtractors::Scalar z_displacement(2); if (apply_dirichlet_bc == true) VectorTools::interpolate_boundary_values(dof_handler_ref, -- 2.39.5