Hello,
I have stripped down my code to an example containing a neo-hookean
solid with displacement boundary conditions. See the attached patch. I
have just tested this with the default configure options. Any
suggestions for improvements are of course welcome.
While porting my code to the libmesh tensor classes, I discovered a
small bug in TypeTensor<T>::assign. See tensor.patch.
Regards
Robert
On Mon, Jan 09, 2012 at 09:02:29AM +1300, Mark Davies wrote:
> Hi Guys,
>
> I'm just getting started with the library and was wondering if there are
> any examples dealing with elastic solids? Thanks in advance for any
> suggestions.
>
> Regards,
> Mark
>
>
> ------------------------------------------------------------------------------
> Ridiculously easy VDI. With Citrix VDI-in-a-Box, you don't need a complex
> infrastructure or vast IT resources to deliver seamless, secure access to
> virtual desktops. With this all-in-one solution, easily deploy virtual
> desktops for less than the cost of PCs and save 60% on VDI infrastructure
> costs. Try it free! http://p.sf.net/sfu/Citrix-VDIinabox
> _______________________________________________
> Libmesh-users mailing list
> libmesh-us...@lists.sourceforge.net
> https://lists.sourceforge.net/lists/listinfo/libmesh-users
Index: examples/fem_system/solid_ex2/solid.in
===================================================================
--- examples/fem_system/solid_ex2/solid.in (Revision 0)
+++ examples/fem_system/solid_ex2/solid.in (Revision 0)
@@ -0,0 +1,21 @@
+schedule/dt = 0.2
+solver/quiet = false
+solver/nonlinear/max_nonlinear_iterations = 100
+solver/nonlinear/relative_step_tolerance = 1.e-3
+solver/nonlinear/relative_residual_tolerance = 1.e-8
+solver/nonlinear/absolute_residual_tolerance = 1.e-8
+solver/quiet = false
+solver/nonlinear/require_reduction = false
+max_linear_iterations = 50000
+initial_linear_tolerance = 1.e-3
+assembly/use_symmetry = true
+bc/displacement = '0 0.0 0.0 0.0 5 nan nan -0.75'
+bc/displacement_penalty = 1e7
+material/neohooke/e_modulus = 10000.0
+material/neohooke/nu = 0.3
+mesh/generation/num_nodes = '4 4 4'
+mesh/generation/origin = '0.0 0.0 0.0'
+mesh/generation/size = '1.5 1.5 1.5'
+duration = 1.0
+output/frequency = 1
+results_directory = ./
Index: examples/fem_system/solid_ex2/solid_system.C
===================================================================
--- examples/fem_system/solid_ex2/solid_system.C (Revision 0)
+++ examples/fem_system/solid_ex2/solid_system.C (Revision 0)
@@ -0,0 +1,346 @@
+#include "getpot.h"
+
+#include "boundary_info.h"
+#include "fe_base.h"
+#include "fem_context.h"
+#include "dof_map.h"
+#include "mesh.h"
+#include "quadrature.h"
+
+#include "equation_systems.h"
+#include "transient_system.h"
+#include "steady_solver.h"
+#include "diff_solver.h"
+#include "newton_solver.h"
+
+#include "sparse_matrix.h"
+#include "numeric_vector.h"
+
+#include "nonlinear_neohooke_cc.h"
+#include "solid_system.h"
+
+// Bring in everything from the libMesh namespace
+using namespace libMesh;
+
+SolidSystem::SolidSystem(EquationSystems& es, const std::string& name,
+ const unsigned int number) :
+ FEMSystem(es, name, number) {
+
+ // Add a time solver. We are just looking at a steady state problem here.
+ this->time_solver = AutoPtr<TimeSolver>(new SteadySolver(*this));
+}
+
+void SolidSystem::save_initial_mesh() {
+ System & aux_sys = this->get_equation_systems().get_system("auxiliary");
+
+ const unsigned int dim = this->get_mesh().mesh_dimension();
+
+ // Loop over all nodes and copy the location from the current system to
+ // the auxiliary system.
+ const MeshBase::const_node_iterator nd_end =
+ this->get_mesh().local_nodes_end();
+ for (MeshBase::const_node_iterator nd = this->get_mesh().local_nodes_begin();
+ nd != nd_end; ++nd) {
+ const Node *node = *nd;
+ for (unsigned int d = 0; d < dim; ++d) {
+ unsigned int source_dof = node->dof_number(this->number(), var[d], 0);
+ unsigned int dest_dof = node->dof_number(aux_sys.number(), undefo_var[d],
+ 0);
+ Real value = this->current_local_solution->el(source_dof);
+ aux_sys.current_local_solution->set(dest_dof, value);
+ }
+ }
+}
+
+void SolidSystem::init_data() {
+ const unsigned int dim = this->get_mesh().mesh_dimension();
+
+ // Get the default order of the used elements. Assumption:
+ // Just one type of elements in the mesh.
+ Order order = this->get_mesh().elem(0)->default_order();
+
+ // Add the node positions as primary variables.
+ var[0] = this->add_variable("x", order);
+ var[1] = this->add_variable("y", order);
+ if (dim == 3)
+ var[2] = this->add_variable("z", order);
+ else
+ var[2] = var[1];
+
+ // Add variables for storing the initial mesh to an auxiliary system.
+ System& aux_sys = this->get_equation_systems().get_system("auxiliary");
+ undefo_var[0] = aux_sys.add_variable("undefo_x", order);
+ undefo_var[1] = aux_sys.add_variable("undefo_y", order);
+ undefo_var[2] = aux_sys.add_variable("undefo_z", order);
+
+ // Set the time stepping options
+ this->deltat = args("schedule/dt", 0.2);
+
+ // Do the parent's initialization after variables are defined
+ FEMSystem::init_data();
+
+// // Tell the system to march velocity forward in time, but
+// // leave p as a constraint only
+// this->time_evolving(var[0]);
+// this->time_evolving(var[1]);
+// if (dim == 3)
+// this->time_evolving(var[2]);
+
+ // Tell the system which variables are containing the node positions
+ set_mesh_system(this);
+
+ this->set_mesh_x_var(var[0]);
+ this->set_mesh_y_var(var[1]);
+ if (dim == 3)
+ this->set_mesh_z_var(var[2]);
+
+ // Fill the variables with the position of the nodes
+ this->mesh_position_get();
+
+ System::reinit();
+
+ // Set some options for the DiffSolver
+ DiffSolver &solver = *(this->time_solver->diff_solver().get());
+ solver.quiet = args("solver/quiet", false);
+ solver.max_nonlinear_iterations = args(
+ "solver/nonlinear/max_nonlinear_iterations", 100);
+ solver.relative_step_tolerance = args(
+ "solver/nonlinear/relative_step_tolerance", 1.e-3);
+ solver.relative_residual_tolerance = args(
+ "solver/nonlinear/relative_residual_tolerance", 1.e-8);
+ solver.absolute_residual_tolerance = args(
+ "solver/nonlinear/absolute_residual_tolerance", 1.e-8);
+ solver.verbose = !args("solver/quiet", false);
+
+ ((NewtonSolver&) solver).require_residual_reduction = args(
+ "solver/nonlinear/require_reduction", false);
+
+ // And the linear solver options
+ solver.max_linear_iterations = args("max_linear_iterations", 50000);
+ solver.initial_linear_tolerance = args("initial_linear_tolerance", 1.e-3);
+}
+
+void SolidSystem::update() {
+ System::update();
+ this->mesh_position_set();
+}
+
+void SolidSystem::init_context(DiffContext &context) {
+ FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
+
+ // Pre-request all the data needed
+ c.element_fe_var[var[0]]->get_JxW();
+ c.element_fe_var[var[0]]->get_phi();
+ c.element_fe_var[var[0]]->get_dphi();
+ c.element_fe_var[var[0]]->get_xyz();
+
+ c.side_fe_var[var[0]]->get_JxW();
+ c.side_fe_var[var[0]]->get_phi();
+ c.side_fe_var[var[0]]->get_xyz();
+}
+
+/**
+ * Compute contribution from internal forces in elem_residual and contribution
from
+ * linearized internal forces (stiffness matrix) in elem_jacobian.
+ */
+bool SolidSystem::element_time_derivative(bool request_jacobian,
+ DiffContext &context) {
+ FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
+
+ // First we get some references to cell-specific data that
+ // will be used to assemble the linear system.
+
+ // Element Jacobian * quadrature weights for interior integration
+ const std::vector<Real> &JxW = c.element_fe_var[var[0]]->get_JxW();
+
+ // The gradients of the shape functions at interior
+ // quadrature points.
+ const std::vector<std::vector<RealGradient> >& dphi =
+ c.element_fe_var[var[0]]->get_dphi();
+
+ // Dimension of the mesh
+ const unsigned int dim = this->get_mesh().mesh_dimension();
+
+ // The number of local degrees of freedom in each variable
+ const unsigned int n_u_dofs = c.dof_indices_var[var[0]].size();
+ libmesh_assert(n_u_dofs == c.dof_indices_var[var[1]].size());
+ if (dim == 3) {
+ libmesh_assert(n_u_dofs == c.dof_indices_var[var[2]].size());
+ }
+
+ unsigned int n_qpoints = c.element_qrule->n_points();
+
+ // Some matrices and vectors for storing the results of the constitutive
+ // law
+ DenseMatrix<Real> stiff;
+ DenseVector<Real> res;
+ VectorValue<Gradient> grad_u;
+
+ // Instantiate the constitutive law
+ NonlinearNeoHookeCurrentConfig material(dphi, args);
+
+ // Just calculate jacobian contribution when we need to
+ material.calculate_linearized_stiffness = request_jacobian;
+
+ // Get a reference to the auxiliary system
+ TransientExplicitSystem& aux_system =
this->get_equation_systems().get_system<
+ TransientExplicitSystem>("auxiliary");
+ std::vector<unsigned int> undefo_index;
+
+ // Assume symmetry of local stiffness matrices
+ bool use_symmetry = args("assembly/use_symmetry", false);
+
+ // Now we will build the element Jacobian and residual.
+ // This must be calculated at each quadrature point by
+ // summing the solution degree-of-freedom values by
+ // the appropriate weight functions.
+ // This class just takes care of the assembly. The matrix of
+ // the jacobian and the residual vector are provided by the
+ // constitutive formulation.
+
+ for (unsigned int qp = 0; qp != n_qpoints; qp++) {
+ // Compute the displacement gradient
+ grad_u(0) = grad_u(1) = grad_u(2) = 0;
+ for (unsigned int d = 0; d < dim; ++d) {
+ std::vector<Number> u_undefo;
+ aux_system.get_dof_map().dof_indices(c.elem, undefo_index,
undefo_var[d]);
+ aux_system.current_local_solution->get(undefo_index, u_undefo);
+ for (unsigned int l = 0; l != n_u_dofs; l++)
+ grad_u(d).add_scaled(dphi[l][qp], u_undefo[l]); // u_current(l)); // -
+ }
+
+ // initialize the constitutive formulation with the current displacement
+ // gradient
+ material.init_for_qp(grad_u, qp);
+
+ // Aquire, scale and assemble residual and stiffness
+ for (unsigned int i = 0; i < n_u_dofs; i++) {
+ res.resize(dim);
+ material.get_residual(res, i);
+ res.scale(JxW[qp]);
+ for (unsigned int ii = 0; ii < dim; ++ii) {
+ c.elem_subresiduals[ii]->operator ()(i) += res(ii);
+ }
+
+ if (request_jacobian && c.elem_solution_derivative) {
+ libmesh_assert(c.elem_solution_derivative == 1.0);
+ for (unsigned int j = (use_symmetry ? i : 0); j < n_u_dofs; j++) {
+ material.get_linearized_stiffness(stiff, i, j);
+ stiff.scale(JxW[qp]);
+ for (unsigned int ii = 0; ii < dim; ++ii) {
+ for (unsigned int jj = 0; jj < dim; ++jj) {
+ c.elem_subjacobians[ii][jj]->operator ()(i, j) += stiff(ii, jj);
+ if (use_symmetry && i != j) {
+ c.elem_subjacobians[ii][jj]->operator ()(j, i) += stiff(jj,
ii);
+ }
+ }
+ }
+ }
+ }
+ }
+ } // end of the quadrature point qp-loop
+
+ return request_jacobian;
+}
+
+bool SolidSystem::side_time_derivative(bool request_jacobian,
+ DiffContext &context) {
+ FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
+
+ // Apply displacement boundary conditions with penalty method
+
+ // Get the current load step
+ Real ratio = this->get_equation_systems().parameters.get<Real>("progress")
+ + 0.001;
+
+ // The BC are stored in the simulation parameters as array containing
sequences of
+ // four numbers: Id of the side for the displacements and three values
describing the
+ // displacement. E.g.: bc/displacement = '5 nan nan -1.0'. This will move
all nodes of
+ // side 5 about 1.0 units down the z-axis while leaving all other directions
unrestricted
+
+ // Get number of BCs to enforce
+ unsigned int num_bc = args.vector_variable_size("bc/displacement");
+ if (num_bc % 4 != 0) {
+ libMesh::err
+ << "ERROR, Odd number of values in displacement boundary condition.\n"
+ << std::endl;
+ libmesh_error();
+ }
+ num_bc /= 4;
+
+ // Loop over all BCs
+ for (unsigned int nbc = 0; nbc < num_bc; nbc++) {
+ // Get IDs of the side for this BC
+ short int positive_boundary_id = args("bc/displacement", 1, nbc * 4);
+ // Check whether this element is on the current boundary
+ short int boundary_id = this->get_mesh().boundary_info->boundary_id(c.elem,
+ c.side);
+
+ // The current side is not on boundary to be restricted
+ if (boundary_id != positive_boundary_id)
+ continue;
+
+ // Read values from configuration file
+ Point diff_value;
+ for (unsigned int d = 0; d < c.dim; ++d) {
+ diff_value(d) = args("bc/displacement", NAN, nbc * 4 + 1 + d);
+ }
+ // Scale according to current load step
+ diff_value *= ratio;
+
+ Real penalty_number = args("bc/displacement_penalty", 1e7);
+
+ FEBase * fe = c.side_fe_var[var[0]];
+ const std::vector<std::vector<Real> > & phi = fe->get_phi();
+ const std::vector<Real>& JxW = fe->get_JxW();
+ const std::vector<Point>& coords = fe->get_xyz();
+
+ unsigned int n_x_dofs = c.dof_indices_var[this->var[0]].size();
+
+ // get mappings for dofs for auxiliary system for original mesh positions
+ const System & auxsys = this->get_equation_systems().get_system(
+ "auxiliary");
+ const DofMap & auxmap = auxsys.get_dof_map();
+ std::vector<unsigned int> undefo_dofs[3];
+ for (unsigned int d = 0; d < c.dim; ++d) {
+ auxmap.dof_indices(c.elem, undefo_dofs[d], undefo_var[d]);
+ }
+
+ for (unsigned int qp = 0; qp < c.side_qrule->n_points(); ++qp) {
+ // calculate coordinates of qp on undeformed mesh
+ Point orig_point;
+ for (unsigned int i = 0; i < n_x_dofs; ++i) {
+ for (unsigned int d = 0; d < c.dim; ++d) {
+ Real orig_val = auxsys.current_solution(undefo_dofs[d][i]);
+ orig_point(d) += phi[i][qp] * orig_val;
+ }
+ }
+
+ // Calculate displacement to be enforced.
+ Point diff = coords[qp] - orig_point - diff_value;
+
+ // Assemble
+ for (unsigned int i = 0; i < n_x_dofs; ++i) {
+ for (unsigned int d1 = 0; d1 < c.dim; ++d1) {
+ if (isnan(diff(d1)))
+ continue;
+ Real val = JxW[qp] * phi[i][qp] * diff(d1) * penalty_number;
+ c.elem_subresiduals[var[d1]]->operator ()(i) += val;
+ }
+ if (request_jacobian) {
+ for (unsigned int j = 0; j < n_x_dofs; ++j) {
+ for (unsigned int d1 = 0; d1 < c.dim; ++d1) {
+ if (isnan(diff(d1)))
+ continue;
+ Real val = JxW[qp] * phi[i][qp] * phi[j][qp] * penalty_number;
+ c.elem_subjacobians[var[d1]][var[d1]]->operator ()(i, j) += val;
+ }
+ }
+ }
+ }
+ }
+ }
+
+ return request_jacobian;
+}
+
Index: examples/fem_system/solid_ex2/nonlinear_neohooke_cc.C
===================================================================
--- examples/fem_system/solid_ex2/nonlinear_neohooke_cc.C (Revision 0)
+++ examples/fem_system/solid_ex2/nonlinear_neohooke_cc.C (Revision 0)
@@ -0,0 +1,150 @@
+#include "nonlinear_neohooke_cc.h"
+
+/**
+ * Return the inverse of the given TypeTensor. Algorithm taken from the tensor
classes
+ * of Ton van den Boogaard (http://tmku209.ctw.utwente.nl/~ton/tensor.html)
+ */
+template <typename T> TypeTensor<T> inv(const TypeTensor<T> &A ) {
+ double Sub11, Sub12, Sub13;
+ Sub11 = A._coords[4]*A._coords[8] - A._coords[5]*A._coords[7];
+ Sub12 = A._coords[3]*A._coords[8] - A._coords[6]*A._coords[5];
+ Sub13 = A._coords[3]*A._coords[7] - A._coords[6]*A._coords[4];
+ double detA = A._coords[0]*Sub11 - A._coords[1]*Sub12 + A._coords[2]*Sub13;
+ libmesh_assert( std::fabs(detA)>1.e-15 );
+
+ TypeTensor<T> Ainv(A);
+
+ Ainv._coords[0] = Sub11/detA;
+ Ainv._coords[1] =
(-A._coords[1]*A._coords[8]+A._coords[2]*A._coords[7])/detA;
+ Ainv._coords[2] = (
A._coords[1]*A._coords[5]-A._coords[2]*A._coords[4])/detA;
+ Ainv._coords[3] = -Sub12/detA;
+ Ainv._coords[4] = (
A._coords[0]*A._coords[8]-A._coords[2]*A._coords[6])/detA;
+ Ainv._coords[5] =
(-A._coords[0]*A._coords[5]+A._coords[2]*A._coords[3])/detA;
+ Ainv._coords[6] = Sub13/detA;
+ Ainv._coords[7] =
(-A._coords[0]*A._coords[7]+A._coords[1]*A._coords[6])/detA;
+ Ainv._coords[8] = (
A._coords[0]*A._coords[4]-A._coords[1]*A._coords[3])/detA;
+
+ return Ainv;
+}
+
+void NonlinearNeoHookeCurrentConfig::init_for_qp(VectorValue<Gradient> &
grad_u, unsigned int qp) {
+ this->current_qp = qp;
+ F.zero();
+ S.zero();
+
+ {
+ RealTensor invF;
+ invF.zero();
+ for (unsigned int i = 0; i < 3; ++i)
+ for (unsigned int j = 0; j < 3; ++j) {
+ invF(i, j) += grad_u(i)(j);
+ }
+ F.add(inv(invF));
+ }
+
+ if (F.det() < -TOLERANCE) {
+ std::cout << "detF < 0" << std::endl;
+ libmesh_error();
+ }
+
+ if (this->calculate_linearized_stiffness) {
+ this->calculate_tangent();
+ }
+
+ this->calculate_stress();
+}
+
+
+
+void NonlinearNeoHookeCurrentConfig::calculate_tangent() {
+ Real mu = E / (2 * (1 + nu));
+ Real lambda = E * nu / ((1 + nu) * (1 - 2 * nu));
+
+ Real detF = F.det();
+
+ C_mat.resize(6, 6);
+ for (unsigned int i = 0; i < 3; ++i) {
+ for (unsigned int j = 0; j < 3; ++j) {
+ if (i == j) {
+ C_mat(i, j) = 2 * mu + lambda;
+ C_mat(i + 3, j + 3) = mu - 0.5 * lambda * (detF
* detF - 1);
+ } else {
+ C_mat(i, j) = lambda * detF * detF;
+ }
+ }
+ }
+}
+
+
+void NonlinearNeoHookeCurrentConfig::calculate_stress() {
+
+ double mu = E / (2.0 * (1.0 + nu));
+ double lambda = E * nu / ((1 + nu) * (1 - 2 * nu));
+
+ Real detF = F.det();
+ RealTensor Ft = F.transpose();
+
+ RealTensor C = Ft * F;
+ RealTensor b = F * Ft;
+ RealTensor identity;
+ identity(0, 0) = 1.0; identity(1, 1) = 1.0; identity(2, 2) = 1.0;
+ RealTensor invC = inv(C);
+
+ S = 0.5 * lambda * (detF * detF - 1) * invC + mu * (identity - invC);
+
+ tau = (F * S) * Ft;
+ sigma = 1/detF * tau;
+}
+
+void NonlinearNeoHookeCurrentConfig::get_residual(DenseVector<Real> &
residuum, unsigned int & i) {
+ B_L.resize(3, 6);
+ DenseVector<Real> sigma_voigt(6);
+
+ this->build_b_0_mat(i, B_L);
+
+ tensor_to_voigt(sigma, sigma_voigt);
+
+ B_L.vector_mult(residuum, sigma_voigt);
+}
+
+void NonlinearNeoHookeCurrentConfig::tensor_to_voigt(const RealTensor &tensor,
DenseVector<Real> &vec) {
+ vec(0) = tensor(0, 0);
+ vec(1) = tensor(1, 1);
+ vec(2) = tensor(2, 2);
+ vec(3) = tensor(0, 1);
+ vec(4) = tensor(1, 2);
+ vec(5) = tensor(0, 2);
+
+}
+
+void
NonlinearNeoHookeCurrentConfig::get_linearized_stiffness(DenseMatrix<Real> &
stiffness, unsigned int & i, unsigned int & j) {
+ stiffness.resize(3, 3);
+
+ double G_IK = (sigma * dphi[i][current_qp]) * dphi[j][current_qp];
+ stiffness(0, 0) += G_IK;
+ stiffness(1, 1) += G_IK;
+ stiffness(2, 2) += G_IK;
+
+ B_L.resize(3, 6);
+ this->build_b_0_mat(i, B_L);
+ B_K.resize(3, 6);
+ this->build_b_0_mat(j, B_K);
+
+ B_L.right_multiply(C_mat);
+ B_L.right_multiply_transpose(B_K);
+ B_L *= 1/F.det();
+
+ stiffness += B_L;
+}
+
+void NonlinearNeoHookeCurrentConfig::build_b_0_mat(int i, DenseMatrix<Real>&
b_0_mat) {
+ for (unsigned int ii = 0; ii < 3; ++ii) {
+ b_0_mat(ii, ii) = dphi[i][current_qp](ii);
+ }
+ b_0_mat(0, 3) = dphi[i][current_qp](1);
+ b_0_mat(1, 3) = dphi[i][current_qp](0);
+ b_0_mat(1, 4) = dphi[i][current_qp](2);
+ b_0_mat(2, 4) = dphi[i][current_qp](1);
+ b_0_mat(0, 5) = dphi[i][current_qp](2);
+ b_0_mat(2, 5) = dphi[i][current_qp](0);
+}
Index: examples/fem_system/solid_ex2/solid_system.h
===================================================================
--- examples/fem_system/solid_ex2/solid_system.h (Revision 0)
+++ examples/fem_system/solid_ex2/solid_system.h (Revision 0)
@@ -0,0 +1,55 @@
+#ifndef SOLID_SYSTEM_H_
+#define SOLID_SYSTEM_H_
+
+#include "fem_system.h"
+#include "auto_ptr.h"
+
+using namespace libMesh;
+
+class SolidSystem: public FEMSystem {
+public:
+ // Constructor
+ SolidSystem(EquationSystems& es, const std::string& name,
+ const unsigned int number);
+
+ // System initialization
+ virtual void init_data();
+
+ // Context initialization
+ virtual void init_context(DiffContext &context);
+
+ // Element residual and jacobian calculations
+ virtual bool element_time_derivative(bool request_jacobian,
+ DiffContext& context);
+
+ // Contributions for adding boundary conditions
+ virtual bool side_time_derivative(bool request_jacobian,
+ DiffContext& context);
+
+ virtual bool eulerian_residual(bool, DiffContext &) {
+ return false;
+ }
+
+ // Simulation parameters
+ GetPot args;
+
+ // Custom Identifier
+ virtual std::string system_type() const {
+ return "Solid";
+ }
+
+ // override method to update mesh also
+ virtual void update();
+
+ // save the undeformed mesh to an auxiliary system
+ void save_initial_mesh();
+
+ // variable numbers of primary variables in the current system
+ unsigned int var[3];
+
+ // variable numbers of primary variables in auxiliary system (for accessing
the
+ // undeformed configuration)
+ unsigned int undefo_var[3];
+};
+
+#endif /* SOLID_SYSTEM_H_ */
Index: examples/fem_system/solid_ex2/nonlinear_neohooke_cc.h
===================================================================
--- examples/fem_system/solid_ex2/nonlinear_neohooke_cc.h (Revision 0)
+++ examples/fem_system/solid_ex2/nonlinear_neohooke_cc.h (Revision 0)
@@ -0,0 +1,65 @@
+#ifndef NONLINEAR_NEOHOOKE_CC_H_
+#define NONLINEAR_NEOHOOKE_CC_H_
+
+#include <dense_vector.h>
+#include <dense_matrix.h>
+#include <vector_value.h>
+#include <tensor_value.h>
+#include <getpot.h>
+
+/**
+ * This class implements a constitutive formulation for an Neo-Hookean elastic
solid
+ * in terms of the current configuration. This implementation is suitable for
computing
+ * large deformations. See e.g. "Nonlinear finite element methods" (P.
Wriggers, 2008,
+ * Springer) for details.
+ */
+class NonlinearNeoHookeCurrentConfig {
+public:
+ NonlinearNeoHookeCurrentConfig(
+ const std::vector<std::vector<RealGradient> >& dphi, GetPot& args) :
+ dphi(dphi) {
+ E = args("material/neohooke/e_modulus", 10000.0);
+ nu = args("material/neohooke/nu", 0.3);
+ }
+
+ /**
+ * Initialize the class for the given displacement gradient at the
+ * specified quadrature point.
+ */
+ void init_for_qp(VectorValue<Gradient> & grad_u, unsigned int qp);
+
+ /**
+ * Return the residual vector for the current state.
+ */
+ void get_residual(DenseVector<Real> & residuum, unsigned int & i);
+
+ /**
+ * Return the stiffness matrix for the current state.
+ */
+ void get_linearized_stiffness(DenseMatrix<Real> & stiffness,
+ unsigned int & i, unsigned int & j);
+
+ /**
+ * Flag to indicate if it is necessary to calculate values for stiffness
+ * matrix during initialization.
+ */
+ bool calculate_linearized_stiffness;
+private:
+ void build_b_0_mat(int i, DenseMatrix<Real>& b_l_mat);
+ void calculate_stress();
+ void calculate_tangent();
+ static void tensor_to_voigt(const RealTensor &tensor, DenseVector<Real>
&vec);
+
+ unsigned int current_qp;
+ const std::vector<std::vector<RealGradient> >& dphi;
+
+ DenseMatrix<Real> C_mat;
+ Real E;
+ Real nu;
+ RealTensor F, S, tau, sigma;
+ DenseMatrix<Real> B_L;
+ DenseMatrix<Real> B_K;
+};
+
+#endif /* NONLINEAR_NEOHOOKE_CC_H_ */
+
Index: examples/fem_system/solid_ex2/solid_ex2.C
===================================================================
--- examples/fem_system/solid_ex2/solid_ex2.C (Revision 0)
+++ examples/fem_system/solid_ex2/solid_ex2.C (Revision 0)
@@ -0,0 +1,137 @@
+#include <string>
+#include <sstream>
+#include <fstream>
+#include <iostream>
+
+#include "libmesh.h"
+#include "mesh.h"
+#include "libmesh_logging.h"
+#include "getpot.h"
+#include <unistd.h>
+#include <stdio.h>
+#include <time.h>
+#include <mesh_generation.h>
+#include <string_to_enum.h>
+#include <transient_system.h>
+#include <vtk_io.h>
+#include <numeric_vector.h>
+#include <time_solver.h>
+#include <equation_systems.h>
+
+using namespace libMesh;
+
+#include "solid_system.h"
+
+void setup(EquationSystems& systems, Mesh& mesh, GetPot& args) {
+
+ const unsigned int dim = mesh.mesh_dimension();
+ libmesh_assert (dim == 2 || dim == 3);
+
+ // Generating Mesh
+ ElemType eltype =
Utility::string_to_enum<ElemType>(args("mesh/generation/element_type", "hex8"));
+ int nx = args("mesh/generation/num_nodes", 4, 0);
+ int ny = args("mesh/generation/num_nodes", 4, 1);
+ int nz = args("mesh/generation/num_nodes", 4, 2);
+ double origx = args("mesh/generation/origin", -1.0, 0);
+ double origy = args("mesh/generation/origin", -1.0, 1);
+ double origz = args("mesh/generation/origin", 0.0, 2);
+ double sizex = args("mesh/generation/size", 2.0, 0);
+ double sizey = args("mesh/generation/size", 2.0, 1);
+ double sizez = args("mesh/generation/size", 2.0, 2);
+ MeshTools::Generation::build_cube(mesh, nx, ny, nz,
+ origx, origx+sizex, origy, origy+sizey, origz, origz+sizez, eltype);
+
+ // Creating Systems
+ SolidSystem& imms = systems.add_system<SolidSystem> ("solid");
+ imms.args = args;
+
+ // Build up auxiliary system
+ ExplicitSystem& aux_sys =
systems.add_system<TransientExplicitSystem>("auxiliary");
+
+ // Initialize the system
+ systems.parameters.set<unsigned int>("phase") = 0;
+ systems.init();
+
+ imms.save_initial_mesh();
+
+ // Fill global solution vector from local ones
+ aux_sys.reinit();
+}
+
+void run_timestepping(EquationSystems& systems, GetPot& args) {
+
+ TransientExplicitSystem& aux_system =
systems.get_system<TransientExplicitSystem>("auxiliary");
+
+ SolidSystem& solid_system = systems.get_system<SolidSystem>("solid");
+
+ AutoPtr<VTKIO> io = AutoPtr<VTKIO>(new VTKIO(systems.get_mesh()));
+
+ Real duration = args("duration", 1.0);
+
+ for (unsigned int t_step = 0; t_step < duration/solid_system.deltat;
t_step++) {
+ // Progress in current phase [0..1]
+ Real progress = t_step * solid_system.deltat / duration;
+ systems.parameters.set<Real>("progress") = progress;
+ systems.parameters.set<unsigned int>("step") = t_step;
+
+ // Update message
+
+ out << "===== Time Step " << std::setw(4) << t_step;
+ out << " (" << std::fixed << std::setprecision(2) << std::setw(6) <<
(progress*100.) << "%)";
+ out << ", time = " << std::setw(7) << solid_system.time;
+ out << " =====" << std::endl;
+
+ // Advance timestep in auxiliary system
+ aux_system.current_local_solution->close();
+ aux_system.old_local_solution->close();
+ *aux_system.older_local_solution = *aux_system.old_local_solution;
+ *aux_system.old_local_solution = *aux_system.current_local_solution;
+
+ out << "Solving Solid" << std::endl;
+ solid_system.solve();
+ aux_system.reinit();
+
+ // Carry out the adaptive mesh refinement/coarsening
+ out << "Doing a reinit of the equation systems" << std::endl;
+ systems.reinit();
+
+ if (t_step % args("output/frequency", 1) == 0) {
+ std::string result;
+ std::stringstream file_name;
+ file_name << args("results_directory", "./") << "fem_";
+ file_name << std::setw(6) << std::setfill('0') << t_step;
+ file_name << ".pvtu";
+
+
+ io->write_equation_systems(file_name.str(), systems);
+ }
+ // Advance to the next timestep in a transient problem
+ out << "Advancing to next step" << std::endl;
+ solid_system.time_solver->advance_timestep();
+ }
+}
+
+
+
+int main(int argc, char** argv) {
+
+ // Initialize libMesh and any dependent libraries
+ LibMeshInit init(argc, argv);
+
+ // read simulation parameters from file
+ GetPot args = GetPot("solid.in");
+
+ // Create System and Mesh
+ Mesh mesh(3);
+ EquationSystems systems(mesh);
+
+ // Create and set systems up
+ setup(systems, mesh, args);
+
+ // run the systems
+ run_timestepping(systems, args);
+
+ out << "Finished calculations" << std::endl;
+ return 0;
+}
+
Index: examples/fem_system/solid_ex2/Makefile
===================================================================
--- examples/fem_system/solid_ex2/Makefile (Revision 0)
+++ examples/fem_system/solid_ex2/Makefile (Revision 0)
@@ -0,0 +1,77 @@
+# $Id: Makefile 5000 2011-12-05 04:49:39Z knezed01 $
+
+
+# The location of the mesh library
+LIBMESH_DIR ?= ../../..
+
+# include the library options determined by configure. This will
+# set the variables INCLUDE and LIBS that we will need to build and
+# link with the library.
+include $(LIBMESH_DIR)/Make.common
+
+
+###############################################################################
+# File management. This is where the source, header, and object files are
+# defined
+
+#
+# source files
+srcfiles := $(wildcard *.C)
+
+#
+# object files
+objects := $(patsubst %.C, %.$(obj-suffix), $(srcfiles))
+###############################################################################
+
+
+
+.PHONY: clean clobber distclean
+
+###############################################################################
+# Target:
+#
+target := ./solid_ex1-$(METHOD)
+
+
+all:: $(target)
+
+# Production rules: how to make the target - depends on library configuration
+$(target): $(objects)
+ @echo "Linking "$@"..."
+ @$(libmesh_CXX) $(libmesh_CPPFLAGS) $(libmesh_CXXFLAGS) $(objects) -o
$@ $(libmesh_LIBS) $(libmesh_LDFLAGS)
+
+
+# Useful rules.
+clean:
+ @rm -f $(objects) *.*vtu *~
+
+clobber:
+ @$(MAKE) clean
+ @rm -f $(target) fem*.*vtu
+
+distclean:
+ @$(MAKE) clobber
+ @rm -f *.o *.g.o *.pg.o
+
+run: $(target)
+ @echo "***************************************************************"
+ @echo "* Running Example " $(LIBMESH_RUN) $(target) $(LIBMESH_OPTIONS)
+ @echo "***************************************************************"
+ @echo " "
+ @$(LIBMESH_RUN) $(target) $(LIBMESH_OPTIONS)
+ @echo " "
+ @echo "***************************************************************"
+ @echo "* Done Running Example " $(target)
+ @echo "***************************************************************"
+
+# include the dependency list
+include .depend
+
+
+#
+# Dependencies
+#
+.depend: $(srcfiles) $(LIBMESH_DIR)/include/*/*.h
+ @$(perl) $(LIBMESH_DIR)/contrib/bin/make_dependencies.pl -I. $(foreach
i, $(wildcard $(LIBMESH_DIR)/include/*), -I$(i)) "-S\$$(obj-suffix)"
$(srcfiles) > .depend
+
+###############################################################################
Index: include/numerics/type_tensor.h
===================================================================
--- include/numerics/type_tensor.h (Revision 5026)
+++ include/numerics/type_tensor.h (Arbeitskopie)
@@ -408,7 +408,7 @@
inline
void TypeTensor<T>::assign (const TypeTensor<T2> &p)
{
- for (unsigned int i=0; i<LIBMESH_DIM; i++)
+ for (unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
_coords[i] = p._coords[i];
}
------------------------------------------------------------------------------
Write once. Port to many.
Get the SDK and tools to simplify cross-platform app development. Create
new or port existing apps to sell to consumers worldwide. Explore the
Intel AppUpSM program developer opportunity. appdeveloper.intel.com/join
http://p.sf.net/sfu/intel-appdev
_______________________________________________
Libmesh-devel mailing list
Libmesh-devel@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/libmesh-devel
