in case the attached file is not present, i copied-pasted the c++ code
below:
************************************** SOURCE CODE STARTS HERE
***********************************
//=================================================================================================
#include <iostream>
#include <algorithm>
#include <cmath>
#include "dense_submatrix.h"
#include "dense_subvector.h"
#include "gmsh_io.h"
#include "vtk_io.h"
#include "xdr_io.h"
#include "plane.h"
#include "boundary_info.h"
#include "libmesh.h"
#include "mesh.h"
#include "mesh_generation.h"
#include "gmv_io.h"
#include "gnuplot_io.h"
#include "linear_implicit_system.h"
#include "equation_systems.h"
#include "fe.h"
#include "quadrature_gauss.h"
#include "dof_map.h"
#include "sparse_matrix.h"
#include "numeric_vector.h"
#include "dense_matrix.h"
#include "dense_vector.h"
#include "perf_log.h"
#include "elem.h"
#include "string_to_enum.h"
#include "getpot.h"
//=================================================================================================
// Global parameters:
// Material parameters:
const Real young =1.0e0,
poisson =0.3;
const Real cc0 =young/(1.0-poisson*poisson),
cc1 =young*poisson/(1.0-poisson*poisson),
cc2 =young/2.0/(1.0+poisson);
// Normal pressure loading:
const Real pressure =1.0;
// Penalty factor for imposing essential BCs:
const Real penalty =1.0e10;
//=================================================================================================
void assemble_navier_cauchy (EquationSystems& es, const std::string&
system_name);
//=================================================================================================
int main (int argc, char* argv[]) {
LibMeshInit init(argc, argv);
// Construct a 2D square mesh:
Mesh mesh(2);
MeshTools::Generation::build_square(mesh, 16, 16, 0.0, 1.0, 0.0, 1.0,
TRI3);
mesh.print_info();
//
PerfLog perf_log("square plate under uniform tension");
// Create a linear implicit system for the Navier-Cauchy PDE of 2D
elastostatics:
EquationSystems system(mesh);
system.add_system<LinearImplicitSystem>("Navier-Cauchy");
// Add displacement components:
system.get_system("Navier-Cauchy").add_variable("u0", FIRST, LAGRANGE);
system.get_system("Navier-Cauchy").add_variable("u1", FIRST, LAGRANGE);
// Attach assembly function:
system.get_system("Navier-Cauchy").attach_assemble_function(assemble_navier_cauchy);
system.init();
system.print_info();
// Solve linear system:
system.get_system("Navier-Cauchy").solve();
// Print the results in Gmsh format:
GmshIO(mesh).write_equation_systems("Kirsch2dElastostatics.msh",
system);
// Exit function normally:
return 0;
}
//=================================================================================================
void assemble_navier_cauchy (EquationSystems& es, const std::string&
system_name) {
libmesh_assert(system_name=="Navier-Cauchy");
PerfLog perf_log("matrix assembly");
const MeshBase& mesh =es.get_mesh();
const unsigned int dim = mesh.mesh_dimension();
LinearImplicitSystem& system
=es.get_system<LinearImplicitSystem>("Navier-Cauchy");
QGauss quad_volu(dim, FIFTH);
QGauss quad_surf(dim-1, FIFTH);
const unsigned int u0_var =system.variable_number("u0");
const unsigned int u1_var =system.variable_number("u1");
FEType fe_u_type =system.variable_type(u0_var);
AutoPtr<FEBase> fe_u_volu(FEBase::build(dim, fe_u_type)); // FEM for
volume integration.
AutoPtr<FEBase> fe_u_surf(FEBase::build(dim, fe_u_type)); // FEM for
surface integration.
fe_u_volu->attach_quadrature_rule(&quad_volu);
fe_u_surf->attach_quadrature_rule(&quad_surf);
const std::vector<Real>& JxW_volu =fe_u_volu->get_JxW();
const std::vector<std::vector<Real> >& phi_u_volu =fe_u_volu->get_phi();
const std::vector<std::vector<RealGradient> >& dphi_u_volu
=fe_u_volu->get_dphi();
const DofMap& dof_map =system.get_dof_map();
DenseMatrix<Number> Ke;
DenseVector<Number> Fe;
DenseSubMatrix<Number> K00(Ke), K01(Ke),
K10(Ke), K11(Ke);
DenseSubVector<Number> F0(Fe), F1(Fe);
std::vector<unsigned int> dof_indices;
std::vector<unsigned int> dof_indices_u0;
std::vector<unsigned int> dof_indices_u1;
for ( MeshBase::const_element_iterator
citer_el=mesh.local_elements_begin();
citer_el!=mesh.local_elements_end(); citer_el++) {
const Elem* elem =*citer_el;
// Initialize corresponding arrays:
perf_log.push("initialization");
dof_map.dof_indices(elem, dof_indices);
dof_map.dof_indices(elem, dof_indices_u0, u0_var);
dof_map.dof_indices(elem, dof_indices_u1, u1_var);
const unsigned int n_dofs =dof_indices.size();
const unsigned int n_u0_dofs =dof_indices_u0.size();
const unsigned int n_u1_dofs =dof_indices_u1.size();
Ke.resize(n_dofs, n_dofs);
Fe.resize(n_dofs);
K00.reposition(u0_var*n_u0_dofs, u0_var*n_u0_dofs, n_u0_dofs,
n_u0_dofs);
K01.reposition(u0_var*n_u0_dofs, u1_var*n_u1_dofs, n_u0_dofs,
n_u1_dofs);
K10.reposition(u1_var*n_u1_dofs, u0_var*n_u0_dofs, n_u1_dofs,
n_u0_dofs);
K11.reposition(u1_var*n_u1_dofs, u1_var*n_u1_dofs, n_u1_dofs,
n_u1_dofs);
F0.reposition(u0_var*n_u0_dofs, n_u0_dofs);
F1.reposition(u1_var*n_u1_dofs, n_u1_dofs);
fe_u_volu->reinit(elem);
perf_log.pop("initialization");
// Stiffness matrix:
perf_log.push("stiffness matrix");
for (unsigned int qp=0; qp<quad_volu.n_points(); qp++) {
for (unsigned int i=0; i<n_u0_dofs; i++)
for (unsigned int j=0; j<n_u0_dofs; j++)
K00(i,j) +=
JxW_volu[qp]*(dphi_u_volu[i][qp](0)*cc0*dphi_u_volu[j][qp](0)+
dphi_u_volu[i][qp](1)*cc2*dphi_u_volu[j][qp](1));
for (unsigned int i=0; i<n_u0_dofs; i++)
for (unsigned int j=0; j<n_u1_dofs; j++)
K01(i,j) +=
JxW_volu[qp]*(dphi_u_volu[i][qp](0)*cc1*dphi_u_volu[j][qp](1)+
dphi_u_volu[i][qp](1)*cc2*dphi_u_volu[j][qp](0));
for (unsigned int i=0; i<n_u1_dofs; i++)
for (unsigned int j=0; j<n_u0_dofs; j++)
K10(i,j) +=
JxW_volu[qp]*(dphi_u_volu[i][qp](1)*cc1*dphi_u_volu[j][qp](0)+
dphi_u_volu[i][qp](0)*cc2*dphi_u_volu[j][qp](1));
for (unsigned int i=0; i<n_u1_dofs; i++)
for (unsigned int j=0; j<n_u1_dofs; j++)
K11(i,j) +=
JxW_volu[qp]*(dphi_u_volu[i][qp](1)*cc0*dphi_u_volu[j][qp](1)+
dphi_u_volu[i][qp](0)*cc2*dphi_u_volu[j][qp](0));
}
perf_log.pop("stiffness matrix");
// Body-force vector:
perf_log.push("rhs vector");
// No-body forces are present, that's is why the following loop is
commented.
/* for (unsigned int qp=0; qp<quad_volu.n_points(); qp++) {
for (unsigned int i=0; i<n_u0_dofs; i++)
F0(i) += JxW_volu[qp]*phi_u_volu[i][qp]*0.0;
for (unsigned int i=0; i<n_u1_dofs; i++)
F1(i) += JxW_volu[qp]*phi_u_volu[i][qp]*0.0;
}*/
perf_log.pop("rhs vector");
// Boundary conditions:
perf_log.push("insert BCs");
for (unsigned int s=0; s<elem->n_sides(); s++) {
if (elem->neighbor(s)!=NULL) continue;
const short int bottom =0, // LIBMESH square plate under uniform
tension in 2D
right =1,
top =2,
left =3;
// Get the boundary ID for side 's'. These are set internally by
build_square():
short int bc_id =mesh.boundary_info->boundary_id(elem,s);
if (bc_id==BoundaryInfo::invalid_id) libmesh_error();
AutoPtr<Elem> side(elem->build_side(s));
const std::vector<Real>& JxW_surf =fe_u_surf->get_JxW();
const std::vector<std::vector<Real> >& phi_u_surf
=fe_u_surf->get_phi();
const std::vector<Point>& eta =fe_u_surf->get_normals();
fe_u_surf->reinit(elem, s);
if (bc_id==left) { // zero u0 Dirichlet BC:
const Real u0_value =0.0;
for (unsigned int ns=0; ns<side->n_nodes(); ns++) {
for (unsigned int ne=0; ne<elem->n_nodes(); ne++) {
if (elem->node(ne) == side->node(ns)) {
K00(ne,ne) += penalty;
F0(ne) += penalty * u0_value;
}
}
}
} else if (bc_id==bottom) { // zero u1 Dirichlet BC:
const Real u1_value =0.0;
for (unsigned int ns=0; ns<side->n_nodes(); ns++) {
for (unsigned int ne=0; ne<elem->n_nodes(); ne++) {
if (elem->node(ne) == side->node(ns)) {
K11(ne,ne) += penalty;
F1(ne) += penalty * u1_value;
}
}
}
} else if (bc_id==right) { // traction (normal pressure) load
Neumann BC:
for (unsigned int qp=0; qp<quad_surf.n_points(); qp++) {
Real p_value =0.0;
for (unsigned int j=0; j<phi_u_surf.size(); j++) {
p_value += phi_u_surf[j][qp] * pressure;
}
for (unsigned int i=0; i<n_u0_dofs; i++)
F0(i) += JxW_surf[qp] * phi_u_surf[i][qp] *
eta[qp](0) * p_value;
for (unsigned int i=0; i<n_u1_dofs; i++)
F1(i) += JxW_surf[qp] * phi_u_surf[i][qp] *
eta[qp](1) * p_value;
}
}
}
perf_log.pop("insert BCs");
// Insert array and vector in general system:
perf_log.push("matrix insertion");
system.matrix->add_matrix(Ke, dof_indices);
system.rhs->add_vector(Fe, dof_indices);
perf_log.pop("matrix insertion");
}
}
//=================================================================================================
************************************** SOURCE CODE ENDS HERE
***********************************
:)
-------------------------------------------------------------------------
This SF.Net email is sponsored by the Moblin Your Move Developer's challenge
Build the coolest Linux based applications with Moblin SDK & win great prizes
Grand prize is a trip for two to an Open Source event anywhere in the world
http://moblin-contest.org/redirect.php?banner_id=100&url=/
_______________________________________________
Libmesh-users mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/libmesh-users
