It seems like the post below didn't go through, so I'm trying again. Apologies if it's a duplicate.
David
---------- Forwarded message ----------
From: David Knezevic <david.kneze...@akselos.com>
Date: Mon, Dec 8, 2014 at 2:46 PM
Subject: Candidate for systems_of_equations_ex7: nonlinear elasticity
To: libmesh-devel@lists.sourceforge.net
Hi all,
There have been several questions on the mailing list about nonlinear
elasticity, so I wanted to provide a new example (systems_of_equations_ex7)
that addresses this. Note that fem_system_ex2 is a neo-Hookean model, but I
think it's quite complicated, so I think a simpler nonlinear elasticity
example would be helpful for users.
I have attached an implementation of a St. Venant-Kirchoff cantilever beam,
i.e. similar to systems_of_equations_ex6 except that it allows for large
deformation (nonlinear strain). This implementation uses
NonlinearImplicitSystem.
The code converges and seems to give reasonable results, but I was hoping
to get a second (or more!) set of eyes on this to check that it's
reasonable. In particular, the convergence seems to be flaky if I increase
forcing_magnitude, which makes me suspect that something is wrong...
Note also that for testing purposes I recommend using a direct solver
("./example-opt -ksp_type preonly -pc_type lu") so that you don't have to
worry about iterative solver convergence (this example is small so LU is
fast enough).
If anyone is available/interested to look over this example in more detail,
please do let me know.
Thanks,
David
/* The libMesh Finite Element Library. */
/* Copyright (C) 2003 Benjamin S. Kirk */
/* This library is free software; you can redistribute it and/or */
/* modify it under the terms of the GNU Lesser General Public */
/* License as published by the Free Software Foundation; either */
/* version 2.1 of the License, or (at your option) any later version. */
/* This library is distributed in the hope that it will be useful, */
/* but WITHOUT ANY WARRANTY; without even the implied warranty of */
/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU */
/* Lesser General Public License for more details. */
/* You should have received a copy of the GNU Lesser General Public */
/* License along with this library; if not, write to the Free Software */
/* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */
// NonlinearImplicitSystem for St-Venant Kirchoff model (large deformation elasticity).
// We consider an elastic material (hence the constitutive law is the same as in
// example 6) but the system is now nonlinear since we don't neglect the nonlinear
// terms in the strain.
//
// Here we consider a loaded cantilever beam. We consider an elastic material
// (hence the constitutive law is the same as in example 6) but we allow large
// deformations, and hence this is a nonlinear model.
//
// We formulate the PDE on the reference geometry (\Omega) as opposed to the
// deformed geometry (\Omega^deformed). As a result (see Ciarlet's 3D elasticity
// book, Theorem 2.6-2) the PDE is given as follows:
//
// \int_\Omega (F Sigma) : \grad v = \int_\Omega f v + \int_\Gamma g v ds
//
// where:
// * F is the deformation gradient, F = I + du/dx (x here refers to reference coordinates).
// * Sigma is the second Piola-Kirchoff stress. We consider St. Venant Kirchoff material here
// so that Sigma_ij = C_ijkl E_kl (as in systems_of_equations_ex6 except now E_kl contains
// a nonlinear term).
// * E_kl is the strain, E_ij = 0.5 * ( u_i,j + u_j,i + u_k,i u_k,j )
// * f = \det(F) f^(deformed).
// * g = det(F) | F^{-T} n| g^(deformed).
//
// In this example we only consider a body load (e.g. gravity), hence we set g = 0.
// C++ include files that we need
#include <iostream>
#include <algorithm>
#include <cmath>
// Various include files needed for the mesh & solver functionality.
#include "libmesh/libmesh.h"
#include "libmesh/mesh.h"
#include "libmesh/mesh_refinement.h"
#include "libmesh/exodusII_io.h"
#include "libmesh/equation_systems.h"
#include "libmesh/fe.h"
#include "libmesh/quadrature_gauss.h"
#include "libmesh/dof_map.h"
#include "libmesh/sparse_matrix.h"
#include "libmesh/numeric_vector.h"
#include "libmesh/dense_matrix.h"
#include "libmesh/dense_vector.h"
#include "libmesh/elem.h"
#include "libmesh/string_to_enum.h"
#include "libmesh/getpot.h"
#include "libmesh/mesh_generation.h"
#include "libmesh/dirichlet_boundaries.h"
#include "libmesh/zero_function.h"
// The nonlinear solver and system we will be using
#include "libmesh/nonlinear_solver.h"
#include "libmesh/nonlinear_implicit_system.h"
#define BOUNDARY_ID_MIN_Z 0
#define BOUNDARY_ID_MIN_Y 1
#define BOUNDARY_ID_MAX_X 2
#define BOUNDARY_ID_MAX_Y 3
#define BOUNDARY_ID_MIN_X 4
#define BOUNDARY_ID_MAX_Z 5
using namespace libMesh;
EquationSystems *_equation_system = NULL;
/**
* Invert a 3x3 matrix. This is needed in the formulation below.
*/
DenseMatrix<Real> invert_3by3_matrix(const DenseMatrix<Real>& original)
{
libmesh_assert(original.m() == 3 && original.n() == 3);
DenseMatrix<Real> inv(3,3);
Real a11 = original(0,0), a12 = original(0,1), a13 = original(0,2);
Real a21 = original(1,0), a22 = original(1,1), a23 = original(1,2);
Real a31 = original(2,0), a32 = original(2,1), a33 = original(2,2);
Real DET = a11*(a33*a22-a32*a23)-a21*(a33*a12-a32*a13)+a31*(a23*a12-a22*a13);
inv(0,0) = a33*a22-a32*a23; inv(0,1) = -(a33*a12-a32*a13); inv(0,2) = a23*a12-a22*a13;
inv(1,0) = -(a33*a21-a31*a23); inv(1,1) = a33*a11-a31*a13; inv(1,2) = -(a23*a11-a21*a13);
inv(2,0) = a32*a21-a31*a22; inv(2,1) = -(a32*a11-a31*a12); inv(2,2) = a22*a11-a21*a12;
inv.scale(1./DET);
return inv;
}
/**
* Kronecker delta function.
*/
Real kronecker_delta(
unsigned int i,
unsigned int j)
{
return i == j ? 1. : 0.;
}
/**
* Evaluate the fourth order tensor (C_ijkl) that relates stress to strain.
*/
Real elasticity_tensor(
Real young_modulus,
Real poisson_ratio,
unsigned int i,
unsigned int j,
unsigned int k,
unsigned int l)
{
// Define the Lame constants
const Real lambda_1 = (young_modulus*poisson_ratio)/((1.+poisson_ratio)*(1.-2.*poisson_ratio));
const Real lambda_2 = young_modulus/(2.*(1.+poisson_ratio));
return lambda_1 * kronecker_delta(i,j) * kronecker_delta(k,l) +
lambda_2 * (kronecker_delta(i,k) * kronecker_delta(j,l) + kronecker_delta(i,l) * kronecker_delta(j,k));
}
/**
* Evaluate the jacobian of the nonlinear system.
*/
void compute_jacobian (const NumericVector<Number>& soln,
SparseMatrix<Number>& jacobian,
NonlinearImplicitSystem& /*sys*/)
{
EquationSystems &es = *_equation_system;
const Real young_modulus = es.parameters.get<Real>("young_modulus");
const Real poisson_ratio = es.parameters.get<Real>("poisson_ratio");
const Real forcing_magnitude = es.parameters.get<Real>("forcing_magnitude");
const MeshBase& mesh = es.get_mesh();
const unsigned int dim = mesh.mesh_dimension();
NonlinearImplicitSystem& system =
es.get_system<NonlinearImplicitSystem>("NonlinearElasticity");
const unsigned int u_var = system.variable_number ("u");
const unsigned int v_var = system.variable_number ("v");
const unsigned int w_var = system.variable_number ("w");
const DofMap& dof_map = system.get_dof_map();
FEType fe_type = dof_map.variable_type(u_var);
AutoPtr<FEBase> fe (FEBase::build(dim, fe_type));
QGauss qrule (dim, FIFTH);
fe->attach_quadrature_rule (&qrule);
AutoPtr<FEBase> fe_face (FEBase::build(dim, fe_type));
QGauss qface(dim-1, FIFTH);
fe_face->attach_quadrature_rule (&qface);
const std::vector<Real>& JxW = fe->get_JxW();
const std::vector<std::vector<Real> >& phi = fe->get_phi();
const std::vector<std::vector<RealGradient> >& dphi = fe->get_dphi();
DenseMatrix<Number> Ke;
std::vector< std::vector< std::unique_ptr< DenseSubMatrix<Number> > > > Ke_var(3);
for(unsigned int var_i=0; var_i<3; var_i++)
{
Ke_var[var_i].resize(3);
for(unsigned int var_j=0; var_j<3; var_j++)
{
std::unique_ptr< DenseSubMatrix<Number> > Ke_submatrix(new DenseSubMatrix<Number>(Ke));
Ke_var[var_i][var_j] = std::move(Ke_submatrix);
}
}
std::vector<unsigned int> dof_indices;
std::vector< std::vector<unsigned int> > dof_indices_var(3);
jacobian.zero();
MeshBase::const_element_iterator el = mesh.active_local_elements_begin();
const MeshBase::const_element_iterator end_el = mesh.active_local_elements_end();
for ( ; el != end_el; ++el)
{
const Elem* elem = *el;
dof_map.dof_indices (elem, dof_indices);
for(unsigned int var=0; var<3; var++)
{
dof_map.dof_indices (elem, dof_indices_var[var], var);
}
const unsigned int n_dofs = dof_indices.size();
const unsigned int n_var_dofs = dof_indices_var[0].size();
fe->reinit (elem);
Ke.resize (n_dofs,n_dofs);
for(unsigned int var_i=0; var_i<3; var_i++)
for(unsigned int var_j=0; var_j<3; var_j++)
{
Ke_var[var_i][var_j]->reposition (var_i*n_var_dofs, var_j*n_var_dofs, n_var_dofs, n_var_dofs);
}
for (unsigned int qp=0; qp<qrule.n_points(); qp++)
{
DenseVector<Number> u_vec(3);
DenseMatrix<Number> grad_u(3,3);
for(unsigned int var_i=0; var_i<3; var_i++)
{
for (unsigned int j=0; j<n_var_dofs; j++)
{
u_vec(var_i) += phi[j][qp]*soln(dof_indices_var[var_i][j]);
}
for(unsigned int var_j=0; var_j<3; var_j++)
{
for (unsigned int j=0; j<n_var_dofs; j++)
{
// Row is variable u1, u2, or u3, column is x, y, or z
grad_u(var_i,var_j) += dphi[j][qp](var_j)*soln(dof_indices_var[var_i][j]);
}
}
}
DenseMatrix<Real> strain_tensor(3,3);
for(unsigned int i=0; i<3; i++)
for(unsigned int j=0; j<3; j++)
{
strain_tensor(i,j) += 0.5 * ( grad_u(i,j) + grad_u(j,i) );
for(unsigned int k=0; k<3; k++)
{
strain_tensor(i,j) += 0.5 * grad_u(k,i)*grad_u(k,j);
}
}
// Define the deformation gradient
DenseMatrix<Number> F(3,3);
F = grad_u;
for(unsigned int var=0; var<3; var++)
{
F(var,var) += 1.;
}
DenseMatrix<Number> stress_tensor(3,3);
for(unsigned int i=0; i<3; i++)
for(unsigned int j=0; j<3; j++)
for(unsigned int k=0; k<3; k++)
for(unsigned int l=0; l<3; l++)
{
stress_tensor(i,j) +=
elasticity_tensor(young_modulus,poisson_ratio,i,j,k,l) *
strain_tensor(k,l);
}
DenseVector<Number> f_vec(3);
f_vec(0) = 0.;
f_vec(1) = 0.;
f_vec(2) = -forcing_magnitude;
for (unsigned int dof_i=0; dof_i<n_var_dofs; dof_i++)
for (unsigned int dof_j=0; dof_j<n_var_dofs; dof_j++)
{
for(unsigned int i=0; i<3; i++)
for(unsigned int j=0; j<3; j++)
for(unsigned int m=0; m<3; m++)
{
(*Ke_var[i][i])(dof_i,dof_j) += JxW[qp] *
( -dphi[dof_j][qp](m) * stress_tensor(m,j) * dphi[dof_i][qp](j) );
}
for(unsigned int i=0; i<3; i++)
for(unsigned int j=0; j<3; j++)
for(unsigned int k=0; k<3; k++)
for(unsigned int l=0; l<3; l++)
{
Real FxC_ijkl = 0.;
for(unsigned int m=0; m<3; m++)
{
FxC_ijkl +=
F(i,m) *
elasticity_tensor(young_modulus,poisson_ratio,m,j,k,l);
}
(*Ke_var[i][k])(dof_i,dof_j) += JxW[qp] *
( -0.5 * FxC_ijkl *
dphi[dof_j][qp](l) *
dphi[dof_i][qp](j)
);
(*Ke_var[i][l])(dof_i,dof_j) += JxW[qp] *
( -0.5 * FxC_ijkl *
dphi[dof_j][qp](k) *
dphi[dof_i][qp](j)
);
for(unsigned int n=0; n<3; n++)
{
(*Ke_var[i][n])(dof_i,dof_j) += JxW[qp] *
( -0.5 * FxC_ijkl *
( dphi[dof_j][qp](k) * grad_u(n,l) +
dphi[dof_j][qp](l) * grad_u(n,k) ) *
dphi[dof_i][qp](j)
);
}
}
// Jacobian terms due to \int_\Omega f_i det(F) v_i dx
// Use the formula for differentiating a determinant:
// d det(F) / dF_ij = det(F) * (inv(F))_ji.
//
// We want d det(F) / du_i, which, via the chain rule, is:
// d det(F) / du_i = d det(F) / dF_kj * dF_kj / du_i.
// But since F_kj = \delta_kj + u_k,j, we obtain:
// d det(F) / du_i = d det(F) / dF_ij * dF_ij / du_i
// = det(F) * (inv(F))_ji * dF_ij / du_i
// Hence the overall term is:
// f_vec(i) * det(F) * (inv(F))_kj * dF_jk / du_j v_i
for(unsigned int i=0; i<3; i++)
for(unsigned int j=0; j<3; j++)
for(unsigned int k=0; k<3; k++)
{
(*Ke_var[i][j])(dof_i,dof_j) += JxW[qp] *
( F.det() * invert_3by3_matrix(F)(j,k) * dphi[dof_j][qp](k) *
f_vec(i) * phi[dof_i][qp] );
}
}
}
dof_map.constrain_element_matrix (Ke, dof_indices);
jacobian.add_matrix (Ke, dof_indices);
}
}
/**
* Evaluate the residual of the nonlinear system.
*/
void compute_residual (const NumericVector<Number>& soln,
NumericVector<Number>& residual,
NonlinearImplicitSystem& /*sys*/)
{
EquationSystems &es = *_equation_system;
const Real young_modulus = es.parameters.get<Real>("young_modulus");
const Real poisson_ratio = es.parameters.get<Real>("poisson_ratio");
const Real forcing_magnitude = es.parameters.get<Real>("forcing_magnitude");
const MeshBase& mesh = es.get_mesh();
const unsigned int dim = mesh.mesh_dimension();
NonlinearImplicitSystem& system =
es.get_system<NonlinearImplicitSystem>("NonlinearElasticity");
const unsigned int u_var = system.variable_number ("u");
const unsigned int v_var = system.variable_number ("v");
const unsigned int w_var = system.variable_number ("w");
const DofMap& dof_map = system.get_dof_map();
FEType fe_type = dof_map.variable_type(u_var);
AutoPtr<FEBase> fe (FEBase::build(dim, fe_type));
QGauss qrule (dim, FIFTH);
fe->attach_quadrature_rule (&qrule);
AutoPtr<FEBase> fe_face (FEBase::build(dim, fe_type));
QGauss qface(dim-1, FIFTH);
fe_face->attach_quadrature_rule (&qface);
const std::vector<Real>& JxW = fe->get_JxW();
const std::vector<std::vector<Real> >& phi = fe->get_phi();
const std::vector<std::vector<RealGradient> >& dphi = fe->get_dphi();
DenseVector<Number> Re;
std::vector< std::unique_ptr< DenseSubVector<Number> > > Re_var(3);
for(unsigned int var=0; var<3; var++)
{
std::unique_ptr< DenseSubVector<Number> > Re_subvector(new DenseSubVector<Number>(Re));
Re_var[var] = std::move(Re_subvector);
}
std::vector<unsigned int> dof_indices;
std::vector< std::vector<unsigned int> > dof_indices_var(3);
residual.zero();
MeshBase::const_element_iterator el = mesh.active_local_elements_begin();
const MeshBase::const_element_iterator end_el = mesh.active_local_elements_end();
for ( ; el != end_el; ++el)
{
const Elem* elem = *el;
dof_map.dof_indices (elem, dof_indices);
for(unsigned int var=0; var<3; var++)
{
dof_map.dof_indices (elem, dof_indices_var[var], var);
}
const unsigned int n_dofs = dof_indices.size();
const unsigned int n_var_dofs = dof_indices_var[0].size();
fe->reinit (elem);
Re.resize (n_dofs);
for(unsigned int var=0; var<3; var++)
{
Re_var[var]->reposition (var*n_var_dofs, n_var_dofs);
}
for (unsigned int qp=0; qp<qrule.n_points(); qp++)
{
DenseVector<Number> u_vec(3);
DenseMatrix<Number> grad_u(3,3);
for(unsigned int var_i=0; var_i<3; var_i++)
{
for (unsigned int j=0; j<n_var_dofs; j++)
{
u_vec(var_i) += phi[j][qp]*soln(dof_indices_var[var_i][j]);
}
for(unsigned int var_j=0; var_j<3; var_j++)
{
for (unsigned int j=0; j<n_var_dofs; j++)
{
// Row is variable u, v, or w column is x, y, or z
grad_u(var_i,var_j) += dphi[j][qp](var_j)*soln(dof_indices_var[var_i][j]);
}
}
}
DenseMatrix<Real> strain_tensor(3,3);
for(unsigned int i=0; i<3; i++)
for(unsigned int j=0; j<3; j++)
{
strain_tensor(i,j) += 0.5 * ( grad_u(i,j) + grad_u(j,i) );
for(unsigned int k=0; k<3; k++)
{
strain_tensor(i,j) += 0.5 * grad_u(k,i)*grad_u(k,j);
}
}
// Define the deformation gradient
DenseMatrix<Number> F(3,3);
F = grad_u;
for(unsigned int var=0; var<3; var++)
{
F(var,var) += 1.;
}
DenseMatrix<Number> stress_tensor(3,3);
for(unsigned int i=0; i<3; i++)
for(unsigned int j=0; j<3; j++)
for(unsigned int k=0; k<3; k++)
for(unsigned int l=0; l<3; l++)
{
stress_tensor(i,j) +=
elasticity_tensor(young_modulus,poisson_ratio,i,j,k,l) *
strain_tensor(k,l);
}
DenseVector<Number> f_vec(3);
f_vec(0) = 0.;
f_vec(1) = 0.;
f_vec(2) = -forcing_magnitude;
for (unsigned int dof_i=0; dof_i<n_var_dofs; dof_i++)
{
for(unsigned int i=0; i<3; i++)
{
for(unsigned int j=0; j<3; j++)
{
Real FxStress_ij = 0.;
for(unsigned int m=0; m<3; m++)
{
FxStress_ij += F(i,m) * stress_tensor(m,j);
}
(*Re_var[i])(dof_i) += JxW[qp] *
( -FxStress_ij * dphi[dof_i][qp](j) );
}
(*Re_var[i])(dof_i) += JxW[qp] *
( f_vec(i) * F.det() * phi[dof_i][qp] );
}
}
}
dof_map.constrain_element_vector (Re, dof_indices);
residual.add_vector (Re, dof_indices);
}
}
int main (int argc, char** argv)
{
LibMeshInit init (argc, argv);
GetPot infile("nonlinear_elasticity.in");
const Real x_length = infile("x_length", 0.);
const Real y_length = infile("y_length", 0.);
const Real z_length = infile("z_length", 0.);
const Real n_elem_x = infile("n_elem_x", 0);
const Real n_elem_y = infile("n_elem_y", 0);
const Real n_elem_z = infile("n_elem_z", 0);
const std::string approx_order = infile("approx_order", "FIRST");
const std::string fe_family = infile("fe_family", "LAGRANGE");
const Real young_modulus = infile("Young_modulus", 1.0);
const Real poisson_ratio = infile("poisson_ratio", 0.3);
const Real forcing_magnitude = infile("forcing_magnitude", 0.001);
const Real nonlinear_abs_tol = infile("nonlinear_abs_tol", 1.e-8);
const Real nonlinear_rel_tol = infile("nonlinear_rel_tol", 1.e-8);
const unsigned int nonlinear_max_its = infile("nonlinear_max_its", 50);
const unsigned int n_solves = infile("n_solves", 10);
const Real force_scaling = infile("force_scaling", 5.0);
Mesh mesh(init.comm());
MeshTools::Generation::build_cube(
mesh,
n_elem_x,
n_elem_y,
n_elem_z,
0., x_length,
0., y_length,
0., z_length,
HEX27);
mesh.print_info();
EquationSystems equation_systems (mesh);
_equation_system = &equation_systems;
NonlinearImplicitSystem& system =
equation_systems.add_system<NonlinearImplicitSystem> ("NonlinearElasticity");
equation_systems.parameters.set<Real> ("nonlinear solver absolute residual tolerance") = nonlinear_abs_tol;
equation_systems.parameters.set<Real> ("nonlinear solver relative residual tolerance") = nonlinear_rel_tol;
equation_systems.parameters.set<unsigned int> ("nonlinear solver maximum iterations") = nonlinear_max_its;
unsigned int u_var = system.add_variable("u",
Utility::string_to_enum<Order> (approx_order),
Utility::string_to_enum<FEFamily>(fe_family));
unsigned int v_var = system.add_variable("v",
Utility::string_to_enum<Order> (approx_order),
Utility::string_to_enum<FEFamily>(fe_family));
unsigned int w_var = system.add_variable("w",
Utility::string_to_enum<Order> (approx_order),
Utility::string_to_enum<FEFamily>(fe_family));
system.nonlinear_solver->residual = compute_residual;
system.nonlinear_solver->jacobian = compute_jacobian;
equation_systems.parameters.set<Real>("young_modulus") = young_modulus;
equation_systems.parameters.set<Real>("poisson_ratio") = poisson_ratio;
equation_systems.parameters.set<Real>("forcing_magnitude") = forcing_magnitude;
// Attach Dirichlet boundary conditions
std::set<boundary_id_type> clamped_boundaries;
clamped_boundaries.insert(BOUNDARY_ID_MIN_X);
std::vector<unsigned int> uvw;
uvw.push_back(u_var);
uvw.push_back(v_var);
uvw.push_back(w_var);
ZeroFunction<Number> zero;
system.get_dof_map().add_dirichlet_boundary
(DirichletBoundary (clamped_boundaries, uvw, &zero));
equation_systems.init();
equation_systems.print_info();
// Provide a loop here so that we can do a sequence of solves
// where solve n gives a good starting guess for solve n+1.
// This "continuation" approach is helpful for solving for
// large values of "forcing_magnitude".
// Set n_solves and force_scaling in nonlinear_elasticity.in,
for(unsigned int count=0; count<n_solves; count++)
{
Real previous_forcing_magnitude =
equation_systems.parameters.get<Real>("forcing_magnitude");
equation_systems.parameters.set<Real>("forcing_magnitude") =
previous_forcing_magnitude*force_scaling;
libMesh::out << "Performing solve " << count << ", forcing_magnitude: "
<< equation_systems.parameters.get<Real>("forcing_magnitude") << std::endl;
system.solve();
libMesh::out << "System solved at nonlinear iteration "
<< system.n_nonlinear_iterations()
<< " , final nonlinear residual norm: "
<< system.final_nonlinear_residual()
<< std::endl << std::endl;
std::stringstream filename;
filename << "solution_" << count << ".exo";
ExodusII_IO (mesh).write_equation_systems(
filename.str(),
equation_systems);
}
return 0;
}
nonlinear_elasticity.in
Description: Binary data
------------------------------------------------------------------------------ Download BIRT iHub F-Type - The Free Enterprise-Grade BIRT Server from Actuate! Instantly Supercharge Your Business Reports and Dashboards with Interactivity, Sharing, Native Excel Exports, App Integration & more Get technology previously reserved for billion-dollar corporations, FREE http://pubads.g.doubleclick.net/gampad/clk?id=164703151&iu=/4140/ostg.clktrk
_______________________________________________ Libmesh-devel mailing list Libmesh-devel@lists.sourceforge.net https://lists.sourceforge.net/lists/listinfo/libmesh-devel
