Hello guys,
I wrote this program for an elasticity axisimmetric problem of concentric
pipes. I took the example 4 of system equations as basis and made some
changes. The problem is that the result is completely wrong: even thought I
had imposed internal and external pressures and self weight in the right
hand side my output in displacements is constant all over the pipe. At this
point I can only think that I've messed up with the matrix assembling.
Someone knows how can I print the stiffness matrix once assembled?
I've attached a document with the weak formulation of the problem (2D) I'm
trying to solve - it's the equation 9 without the thermal contribute.
Thanks,
rod12
Here it goes:
// C++ include files that we need
#include <iostream>
#include <algorithm>
#include <math.h>
// libMesh includes
#include "libmesh/libmesh.h"
#include "libmesh/mesh.h"
#include "libmesh/mesh_generation.h"
#include "libmesh/exodusII_io.h"
#include "libmesh/gnuplot_io.h"
#include "libmesh/linear_implicit_system.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_submatrix.h"
#include "libmesh/dense_vector.h"
#include "libmesh/dense_subvector.h"
#include "libmesh/perf_log.h"
#include "libmesh/elem.h"
#include "libmesh/boundary_info.h"
#include "libmesh/zero_function.h"
#include "libmesh/dirichlet_boundaries.h"
#include "libmesh/string_to_enum.h"
#include "libmesh/getpot.h"
//Handling errors
void error (std::string & str){
std::cerr<<str<<std::endl;
exit(1);
}
using namespace libMesh;
Real mu (const Real x);
Real lambda (const Real x);
Real rho (const Real x);
Real PressurexR (const Real x,const Real y);
// Bring in everything from the libMesh namespace
// Matrix and right-hand side assemble
void assemble_elasticity(EquationSystems& es,
const std::string& system_name);
// Begin the main program.
int main (int argc, char** argv)
{
// Initialize libMesh and any dependent libaries
LibMeshInit init (argc, argv);
GetPot command_line (argc, argv);
// Initialize the cantilever mesh
const unsigned int dim = 2;
// Skip this 2D example if libMesh was compiled as 1D-only.
//libmesh_example_assert(dim <= LIBMESH_DIM, "2D support");
int psr = 50;
if ( command_line.search(1, "-nx") )
psr = command_line.next(psr);
int psl = 50;
if ( command_line.search(1, "-ny") )
psl = command_line.next(psl);
const Real R_i = 0.12;
const Real R_e = 0.30;
const Real L = 10.0;
// Create a 2D mesh distributed across the default MPI communicator.
Mesh mesh(init.comm(), dim);
MeshTools::Generation::build_square (mesh,
psr, psl,
R_i, R_e,
0., L,
QUAD9);
// Print information about the mesh to the screen.
mesh.print_info();
// Create an equation systems object.
EquationSystems equation_systems (mesh);
// Declare the system and its variables.
// Create a system named "Elasticity"
LinearImplicitSystem& system =
equation_systems.add_system<LinearImplicitSystem> ("Elasticity");
// Add two displacement variables, u and v, to the system
unsigned int u_var = system.add_variable("u", SECOND, LAGRANGE);
unsigned int v_var = system.add_variable("v", SECOND, LAGRANGE);
system.attach_assemble_function (assemble_elasticity);
// Construct a Dirichlet boundary condition object
// We impose a "clamped" boundary condition on the
// "lowear" and "upper" boundaries, i.e. bc_id = 0,2
std::set<boundary_id_type> boundary_ids;
boundary_ids.insert(0);
boundary_ids.insert(2);
// Create a vector storing the variable numbers which the BC applies to
std::vector<unsigned int> variables(2);
variables[0] = u_var; variables[1] = v_var;
// Create a ZeroFunction to initialize dirichlet_bc
ZeroFunction<> zf;
DirichletBoundary dirichlet_bc(boundary_ids,
variables,
&zf);
// We must add the Dirichlet boundary condition _before_
// we call equation_systems.init()
system.get_dof_map().add_dirichlet_boundary(dirichlet_bc);
// Initialize the data structures for the equation system.
equation_systems.init();
// Print information about the system to the screen.
equation_systems.print_info();
// Solve the system
system.solve();
// Plot the solution
#ifdef LIBMESH_HAVE_EXODUS_API
ExodusII_IO
(mesh).write_equation_systems("displacement.e",equation_systems);
#endif // #ifdef LIBMESH_HAVE_EXODUS_API
// All done.
return 0;
}
Real mu (const Real x){
//steel inner and outer constants
const Real E_pipe_mod = 2e11;
const Real poisson_pipe_mod = 0.3;
//mu pipe
const Real mu_p = E_pipe_mod/(2*(1+poisson_pipe_mod));
//insulation constants
const Real E_ins_mod = 5e9;
const Real poisson_ins_mod = 0.4;
const Real mu_i = E_ins_mod/(2*(1+poisson_ins_mod));
//overall internal and external radius
const Real R_i = 0.12;
const Real R_e = 0.30;
//internal and external steel pipe thickness
const Real ti = 0.020;
const Real te = 0.018;
//insulation thickness
if (x >= R_i && x <= R_i+ti )
return mu_p;
else if (x > R_i + ti && x < R_e-te )
return mu_i;
else if (x >= R_e - te && x <= R_e )
return mu_p;
else{
std::string s1("Function mu called for x out of domain bonds!");
error(s1);
}
}
Real lambda (const Real x){
//steel inner and outer constants
const Real E_pipe_mod = 2e11;
const Real poisson_pipe_mod = 0.3;
//mu pipe
const Real lambda_p =
(E_pipe_mod*poisson_pipe_mod)/((1+poisson_pipe_mod)*(1-2*poisson_pipe_mod));
//insulation constants
const Real E_ins_mod = 5e9;
const Real poisson_ins_mod = 0.4;
const Real lambda_i =
(E_ins_mod*poisson_ins_mod)/((1+poisson_ins_mod)*(1-2*poisson_ins_mod));
//overall internal and external radius
const Real R_i = 0.12;
const Real R_e = 0.30;
//internal and external steel pipe thickness
const Real ti = 0.020;
const Real te = 0.018;
//insulation thickness
if (x >= R_i && x <= R_i+ti )
return lambda_p;
else if (x > R_i + ti && x < R_e-te )
return lambda_i;
else if (x >= R_e - te && x <= R_e )
return lambda_p;
else{
std::string s1("Function lambda called for x out of domain bonds!");
error(s1);
}
}
Real rho (const Real x){
//steel inner and outer constants
const Real rho_pipe = 8e3;
const Real rho_ins = 1.3e3;
//overall internal and external radius
const Real R_i = 0.12;
const Real R_e = 0.30;
//internal and external steel pipe thickness
const Real ti = 0.020;
const Real te = 0.018;
//insulation thickness
if (x >= R_i && x <= R_i+ti )
return rho_pipe;
else if (x > R_i + ti && x < R_e-te )
return rho_ins;
else if (x >= R_e - te && x <= R_e )
return rho_pipe;
else{
std::string s1("Function rho called for x out of domain bonds!");
error(s1);
}
}
Real PressurexR (const Real x,const Real y = 0.){
const Real R_i = 0.12;
const Real R_e = 0.30;
const Real L = 10;
//linear change on external temperature with gradient equal to the max
//verified in ocean conditions
if (x == R_i)
return R_i*(104e6 - 10000*y);
else if (x == R_e)
return R_e*(1.5e7 - 10500*y);
else{
std::cerr<<"x = "<<x<<"y = "<<y<<std::endl;
std::string s1("Function PressurexR called for x out of domain
bonds!");
error(s1);
}
}
void assemble_elasticity(EquationSystems& es,
const std::string& system_name)
{
//making sure we will assemble the right system
libmesh_assert_equal_to (system_name, "Elasticity");
//getting the mesh reference
const MeshBase& mesh = es.get_mesh();
//getting dimension of the system
const unsigned int dim = mesh.mesh_dimension();
LinearImplicitSystem& system =
es.get_system<LinearImplicitSystem>("Elasticity");
//getting the variable numbers
const unsigned int u_var = system.variable_number ("u");
const unsigned int v_var = system.variable_number ("v");
//getting dof map and the fe type (both variables have the same type
otherwise need to be changed)
const DofMap& dof_map = system.get_dof_map();
FEType fe_type = dof_map.variable_type(0);
//having a "dynamic" pointer of the prescribed type, getting default
quadrature order
AutoPtr<FEBase> fe (FEBase::build(dim, fe_type));
QGauss qrule (dim, fe_type.default_quadrature_order());
fe->attach_quadrature_rule (&qrule);
//same for boundary elements
AutoPtr<FEBase> fe_face (FEBase::build(dim, fe_type));
QGauss qface(dim-1, fe_type.default_quadrature_order());
fe_face->attach_quadrature_rule (&qface);
//jacobian and quadrature weights for numeric integration (taken in the
quadrature points)
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();
const std::vector<Point>& q_point = fe->get_xyz();
DenseMatrix<Number> Ke;
DenseVector<Number> Fe;
//Submatrix for handling
DenseSubMatrix<Number>
Kuu(Ke), Kuv(Ke),
Kvu(Ke), Kvv(Ke);
DenseSubVector<Number>
Fu(Fe),
Fv(Fe);
std::vector<dof_id_type> dof_indices;
std::vector<dof_id_type> dof_indices_u;
std::vector<dof_id_type> dof_indices_v;
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);
dof_map.dof_indices (elem, dof_indices_u, u_var);
dof_map.dof_indices (elem, dof_indices_v, v_var);
const unsigned int n_dofs = dof_indices.size();
const unsigned int n_u_dofs = dof_indices_u.size();
const unsigned int n_v_dofs = dof_indices_v.size();
fe->reinit (elem);
Ke.resize (n_dofs, n_dofs);
Fe.resize (n_dofs);
Kuu.reposition (u_var*n_u_dofs, u_var*n_u_dofs, n_u_dofs, n_u_dofs);
Kuv.reposition (u_var*n_u_dofs, v_var*n_u_dofs, n_u_dofs, n_v_dofs);
Kvu.reposition (v_var*n_v_dofs, u_var*n_v_dofs, n_v_dofs, n_u_dofs);
Kvv.reposition (v_var*n_v_dofs, v_var*n_v_dofs, n_v_dofs, n_v_dofs);
Fu.reposition (u_var*n_u_dofs, n_u_dofs);
Fv.reposition (v_var*n_u_dofs, n_v_dofs);
for (unsigned int qp=0; qp<qrule.n_points(); qp++)
{
const Real x = q_point[qp](0);
const Real y = q_point[qp](1);
for (unsigned int i=0; i<n_u_dofs; i++){
for (unsigned int j=0; j<n_u_dofs; j++)
{
//See the weak formulation
Kuu(i,j) += JxW[qp]*2*mu(x)*x*dphi[j][qp](0)*dphi[i][qp](0);
Kuu(i,j) += JxW[qp]*2*mu(x)*phi[j][qp]*phi[i][qp]/x;
Kuu(i,j) += JxW[qp]*2*mu(x)*x*dphi[j][qp](1)*dphi[i][qp](1);
Kuu(i,j) += JxW[qp]*lambda(x)*x*dphi[j][qp](0)*dphi[i][qp](0);
Kuu(i,j) += JxW[qp]*lambda(x)*dphi[j][qp](0)*phi[i][qp];
Kuu(i,j) += JxW[qp]*lambda(x)*phi[j][qp]*dphi[i][qp](0);
Kuu(i,j) += JxW[qp]*lambda(x)*phi[j][qp]*phi[i][qp]/x;
}
}
for (unsigned int i=0; i<n_u_dofs; i++)
for (unsigned int j=0; j<n_v_dofs; j++)
{
Kuv(i,j) += JxW[qp]*2*mu(x)*x*dphi[j][qp](0)*dphi[i][qp](1);
Kuv(i,j) += JxW[qp]*lambda(x)*x*dphi[j][qp](1)*dphi[i][qp](0);
Kuv(i,j) += JxW[qp]*lambda(x)*dphi[j][qp](1)*phi[i][qp];
}
for (unsigned int i=0; i<n_v_dofs; i++)
for (unsigned int j=0; j<n_u_dofs; j++)
{
Kvu(i,j) += JxW[qp]*2*mu(x)*x*dphi[j][qp](1)*dphi[i][qp](0);
Kvu(i,j) +=
JxW[qp]*lambda(x)*x*dphi[j][qp](0)*dphi[i][qp](1);
Kvu(i,j) += JxW[qp]*lambda(x)*phi[j][qp]*dphi[j][qp](1);
}
for (unsigned int i=0; i<n_v_dofs; i++)
for (unsigned int j=0; j<n_v_dofs; j++)
{
Kvv(i,j) +=
JxW[qp]*(2*mu(x)+lambda(x))*x*dphi[j][qp](1)*dphi[j][qp](1);
Kvv(i,j) += JxW[qp]*2*mu(x)*x*dphi[j][qp](0)*dphi[j][qp](0);
}
//body load
for (unsigned int i=0; i<n_v_dofs; i++)
Fv(i) += -JxW[qp]*10*x*(1000-rho(x))*phi[i][qp];
}
for (unsigned int side=0; side<elem->n_sides(); side++)
if (elem->neighbor(side) == NULL)
{
const std::vector<std::vector<Real> >& phi_face =
fe_face->get_phi();
const std::vector<Real>& JxW_face = fe_face->get_JxW();
const std::vector<Point>& q_point_face = fe_face->get_xyz();
fe_face->reinit(elem, side);
// Apply pressure on lateral sides (Newmann)
if( !mesh.boundary_info->has_boundary_id (elem, side, 0) &&
!mesh.boundary_info->has_boundary_id (elem, side, 2) )
{
for (unsigned int qp=0; qp<qface.n_points(); qp++)
{
const Real x = q_point_face[qp](0);
const Real y = q_point_face[qp](1);
for (unsigned int i=0; i<n_u_dofs; i++)
{
const Real R_i = 0.12;
const Real R_e = 0.30;
if (mesh.boundary_info->has_boundary_id (elem, side,
1) || mesh.boundary_info->has_boundary_id (elem, side, 3)){
std::cout<<"(x,y)=("<<x<<","<<y<<")"<<std::endl;
std::cout<<"PxR=("<<PressurexR(x,y)<<")"<<std::endl;
Fu(i) +=
-JxW_face[qp]*PressurexR(x,y)*phi_face[i][qp];
}
}
}
}
}
dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices);
system.matrix->add_matrix (Ke, dof_indices);
system.rhs->add_vector (Fe, dof_indices);
}
}
------------------------------------------------------------------------------
Slashdot TV.
Video for Nerds. Stuff that matters.
http://tv.slashdot.org/
_______________________________________________
Libmesh-users mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/libmesh-users
