The source file nonlinear_neohooke_cc.h with comments:
#ifndef NONLINEAR_NEOHOOKE_CC_H_
#define NONLINEAR_NEOHOOKE_CC_H_
#include "libmesh/dense_vector.h"
#include "libmesh/dense_matrix.h"
#include "libmesh/vector_value.h"
#include "libmesh/tensor_value.h"
#include "libmesh/getpot.h"
using namespace libMesh;
/**
* 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_in, GetPot& args) :
dphi(dphi_in) {
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_ */
The source file solid_system.h with comments:
#ifndef SOLID_SYSTEM_H_
#define SOLID_SYSTEM_H_
#include "libmesh/fem_system.h"
#include "libmesh/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_ */
The source file fem_system_ex2.C with comments:
#include "libmesh/equation_systems.h"
#include "libmesh/getpot.h"
#include "libmesh/libmesh.h"
#include "libmesh/libmesh_logging.h"
#include "libmesh/mesh.h"
#include "libmesh/mesh_generation.h"
#include "libmesh/numeric_vector.h"
#include "libmesh/string_to_enum.h"
#include "libmesh/time_solver.h"
#include "libmesh/transient_system.h"
#include "libmesh/vtk_io.h"
#include <cstdio>
#include <ctime>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <sstream>
#include <string>
#include <unistd.h>
using namespace libMesh;
#include "solid_system.h"
void setup(EquationSystems& systems, Mesh& mesh, GetPot& args)
{
const unsigned int dim = mesh.mesh_dimension();
We currently invert tensors with the assumption that they're 3x3
libmesh_assert (dim == 3);
Generating Mesh
ElemType eltype = Utility::string_to_enum<ElemType>(args("mesh/generation/element_type", "hex8"));
int nx = args("mesh/generation/num_elem", 4, 0);
int ny = args("mesh/generation/num_elem", 4, 1);
int nz = dim > 2 ? args("mesh/generation/num_elem", 4, 2) : 0;
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)
{
Skip this example if we do not meet certain requirements
#ifndef LIBMESH_HAVE_VTK
libmesh_example_assert(false, "--enable-vtk");
#endif
Initialize libMesh and any dependent libraries
LibMeshInit init(argc, argv);
Threaded assembly doesn't currently work with the moving mesh
code.
We'll skip this example for now.
if (libMesh::n_threads() > 1)
{
std::cout << "We skip fem_system_ex2 when using threads.\n"
<< std::endl;
return 0;
}
read simulation parameters from file
GetPot args = GetPot("solid.in");
Create System and Mesh
int dim = args("mesh/generation/dimension", 3);
libmesh_example_assert(dim <= LIBMESH_DIM, "3D support");
Mesh mesh(dim);
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;
}
The source file nonlinear_neohooke_cc.C with comments:
#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) {
#if LIBMESH_USE_COMPLEX_NUMBERS
invF(i, j) += grad_u(i)(j).real();
#else
invF(i, j) += grad_u(i)(j);
#endif
}
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);
}
The source file solid_system.C with comments:
#include "libmesh/boundary_info.h"
#include "libmesh/diff_solver.h"
#include "libmesh/dof_map.h"
#include "libmesh/equation_systems.h"
#include "libmesh/fe_base.h"
#include "libmesh/fem_context.h"
#include "libmesh/getpot.h"
#include "libmesh/mesh.h"
#include "libmesh/newton_solver.h"
#include "libmesh/numeric_vector.h"
#include "libmesh/quadrature.h"
#include "libmesh/sparse_matrix.h"
#include "libmesh/steady_solver.h"
#include "libmesh/transient_system.h"
#include "nonlinear_neohooke_cc.h"
#include "solid_system.h"
Solaris Studio has no NAN
#ifdef __SUNPRO_CC
#define NAN (1.0/0.0)
#endif
Bring in everything from the libMesh namespace
using namespace libMesh;
SolidSystem::SolidSystem(EquationSystems& es, const std::string& name_in,
const unsigned int number_in) :
FEMSystem(es, name_in, number_in) {
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);
Number 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().elements_begin()))->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
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
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<dof_id_type> 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
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
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);
The current side may not be on the boundary to be restricted
if (!this->get_mesh().boundary_info->has_boundary_id
(c.elem,c.side,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<dof_id_type> 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) {
Number orig_val = auxsys.current_solution(undefo_dofs[d][i]);
#if LIBMESH_USE_COMPLEX_NUMBERS
orig_point(d) += phi[i][qp] * orig_val.real();
#else
orig_point(d) += phi[i][qp] * orig_val;
#endif
}
}
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 (libmesh_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 (libmesh_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;
}
The source file nonlinear_neohooke_cc.h without comments:
#ifndef NONLINEAR_NEOHOOKE_CC_H_
#define NONLINEAR_NEOHOOKE_CC_H_
#include "libmesh/dense_vector.h"
#include "libmesh/dense_matrix.h"
#include "libmesh/vector_value.h"
#include "libmesh/tensor_value.h"
#include "libmesh/getpot.h"
using namespace libMesh;
/**
* 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_in, GetPot& args) :
dphi(dphi_in) {
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_ */
The source file solid_system.h without comments:
#ifndef SOLID_SYSTEM_H_
#define SOLID_SYSTEM_H_
#include "libmesh/fem_system.h"
#include "libmesh/auto_ptr.h"
using namespace libMesh;
class SolidSystem: public FEMSystem {
public:
SolidSystem(EquationSystems& es, const std::string& name,
const unsigned int number);
virtual void init_data();
virtual void init_context(DiffContext &context);
virtual bool element_time_derivative(bool request_jacobian,
DiffContext& context);
virtual bool side_time_derivative(bool request_jacobian,
DiffContext& context);
virtual bool eulerian_residual(bool, DiffContext &) {
return false;
}
GetPot args;
virtual std::string system_type() const {
return "Solid";
}
virtual void update();
void save_initial_mesh();
unsigned int var[3];
unsigned int undefo_var[3];
};
#endif /* SOLID_SYSTEM_H_ */
The source file fem_system_ex2.C without comments:
#include "libmesh/equation_systems.h"
#include "libmesh/getpot.h"
#include "libmesh/libmesh.h"
#include "libmesh/libmesh_logging.h"
#include "libmesh/mesh.h"
#include "libmesh/mesh_generation.h"
#include "libmesh/numeric_vector.h"
#include "libmesh/string_to_enum.h"
#include "libmesh/time_solver.h"
#include "libmesh/transient_system.h"
#include "libmesh/vtk_io.h"
#include <cstdio>
#include <ctime>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <sstream>
#include <string>
#include <unistd.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 == 3);
ElemType eltype = Utility::string_to_enum<ElemType>(args("mesh/generation/element_type", "hex8"));
int nx = args("mesh/generation/num_elem", 4, 0);
int ny = args("mesh/generation/num_elem", 4, 1);
int nz = dim > 2 ? args("mesh/generation/num_elem", 4, 2) : 0;
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);
SolidSystem& imms = systems.add_system<SolidSystem> ("solid");
imms.args = args;
ExplicitSystem& aux_sys = systems.add_system<TransientExplicitSystem>("auxiliary");
systems.parameters.set<unsigned int>("phase") = 0;
systems.init();
imms.save_initial_mesh();
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++) {
Real progress = t_step * solid_system.deltat / duration;
systems.parameters.set<Real>("progress") = progress;
systems.parameters.set<unsigned int>("step") = t_step;
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;
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();
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);
}
out << "Advancing to next step" << std::endl;
solid_system.time_solver->advance_timestep();
}
}
int main(int argc, char** argv)
{
#ifndef LIBMESH_HAVE_VTK
libmesh_example_assert(false, "--enable-vtk");
#endif
LibMeshInit init(argc, argv);
if (libMesh::n_threads() > 1)
{
std::cout << "We skip fem_system_ex2 when using threads.\n"
<< std::endl;
return 0;
}
GetPot args = GetPot("solid.in");
int dim = args("mesh/generation/dimension", 3);
libmesh_example_assert(dim <= LIBMESH_DIM, "3D support");
Mesh mesh(dim);
EquationSystems systems(mesh);
setup(systems, mesh, args);
run_timestepping(systems, args);
out << "Finished calculations" << std::endl;
return 0;
}
The source file nonlinear_neohooke_cc.C without comments:
#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) {
#if LIBMESH_USE_COMPLEX_NUMBERS
invF(i, j) += grad_u(i)(j).real();
#else
invF(i, j) += grad_u(i)(j);
#endif
}
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);
}
The source file solid_system.C without comments:
#include "libmesh/boundary_info.h"
#include "libmesh/diff_solver.h"
#include "libmesh/dof_map.h"
#include "libmesh/equation_systems.h"
#include "libmesh/fe_base.h"
#include "libmesh/fem_context.h"
#include "libmesh/getpot.h"
#include "libmesh/mesh.h"
#include "libmesh/newton_solver.h"
#include "libmesh/numeric_vector.h"
#include "libmesh/quadrature.h"
#include "libmesh/sparse_matrix.h"
#include "libmesh/steady_solver.h"
#include "libmesh/transient_system.h"
#include "nonlinear_neohooke_cc.h"
#include "solid_system.h"
#ifdef __SUNPRO_CC
#define NAN (1.0/0.0)
#endif
using namespace libMesh;
SolidSystem::SolidSystem(EquationSystems& es, const std::string& name_in,
const unsigned int number_in) :
FEMSystem(es, name_in, number_in) {
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();
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);
Number 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();
Order order = (*(this->get_mesh().elements_begin()))->default_order();
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];
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);
this->deltat = args("schedule/dt", 0.2);
FEMSystem::init_data();
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]);
this->mesh_position_get();
System::reinit();
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);
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);
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);
const std::vector<Real> &JxW = c.element_fe_var[var[0]]->get_JxW();
const std::vector<std::vector<RealGradient> >& dphi =
c.element_fe_var[var[0]]->get_dphi();
const unsigned int dim = this->get_mesh().mesh_dimension();
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();
DenseMatrix<Real> stiff;
DenseVector<Real> res;
VectorValue<Gradient> grad_u;
NonlinearNeoHookeCurrentConfig material(dphi, args);
material.calculate_linearized_stiffness = request_jacobian;
TransientExplicitSystem& aux_system = this->get_equation_systems().get_system<
TransientExplicitSystem>("auxiliary");
std::vector<dof_id_type> undefo_index;
bool use_symmetry = args("assembly/use_symmetry", false);
for (unsigned int qp = 0; qp != n_qpoints; qp++) {
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)); // -
}
material.init_for_qp(grad_u, qp);
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);
Real ratio = this->get_equation_systems().parameters.get<Real>("progress")
+ 0.001;
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;
for (unsigned int nbc = 0; nbc < num_bc; nbc++) {
short int positive_boundary_id = args("bc/displacement", 1, nbc * 4);
if (!this->get_mesh().boundary_info->has_boundary_id
(c.elem,c.side,positive_boundary_id))
continue;
Point diff_value;
for (unsigned int d = 0; d < c.dim; ++d) {
diff_value(d) = args("bc/displacement", NAN, nbc * 4 + 1 + d);
}
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();
const System & auxsys = this->get_equation_systems().get_system(
"auxiliary");
const DofMap & auxmap = auxsys.get_dof_map();
std::vector<dof_id_type> 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) {
Point orig_point;
for (unsigned int i = 0; i < n_x_dofs; ++i) {
for (unsigned int d = 0; d < c.dim; ++d) {
Number orig_val = auxsys.current_solution(undefo_dofs[d][i]);
#if LIBMESH_USE_COMPLEX_NUMBERS
orig_point(d) += phi[i][qp] * orig_val.real();
#else
orig_point(d) += phi[i][qp] * orig_val;
#endif
}
}
Point diff = coords[qp] - orig_point - diff_value;
for (unsigned int i = 0; i < n_x_dofs; ++i) {
for (unsigned int d1 = 0; d1 < c.dim; ++d1) {
if (libmesh_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 (libmesh_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;
}
The console output of the program:
***************************************************************
* Running Example fem_system_ex2:
* mpirun -np 2 example-devel -pc_type bjacobi -sub_pc_type ilu -sub_pc_factor_levels 4 -sub_pc_factor_zeropivot 0 -ksp_right_pc -log_summary
***************************************************************
*** Warning, This code is untested, experimental, or likely to see future API changes: ../../../include/libmesh/diff_physics.h, line 371, compiled Feb 5 2013 at 13:34:45 ***
*** Warning, This code is untested, experimental, or likely to see future API changes: ../../../include/libmesh/vtk_io.h, line 178, compiled Feb 5 2013 at 13:34:44 ***
===== Time Step 0 ( 0.00%), time = 0.00 =====
Solving Solid
Assembling the System
Nonlinear Residual: 3691.41
Linear solve starting, tolerance 0.00
Linear solve finished, step 9, residual 0.00
Trying full Newton step
Current Residual: 0.19
Nonlinear solver current_residual 0.19 > 0.00
Nonlinear solver relative residual 0.00 > 0.00
Nonlinear solver converged, step 0, relative step size 0.00 < 0.00
Doing a reinit of the equation systems
Advancing to next step
===== Time Step 1 ( 20.00%), time = 0.00 =====
Solving Solid
Assembling the System
Nonlinear Residual: 738516.55
Linear solve starting, tolerance 0.00
Linear solve finished, step 9, residual 0.00
Trying full Newton step
Current Residual: 44.89
Nonlinear step: |du|/|u| = 0.06, |du| = 1.04
Assembling the System
Nonlinear Residual: 44.89
Linear solve starting, tolerance 0.00
Linear solve finished, step 9, residual 0.00
Trying full Newton step
Current Residual: 1.52
Nonlinear solver current_residual 1.52 > 0.00
Nonlinear solver relative residual 0.00 > 0.00
Nonlinear solver converged, step 1, relative step size 0.00 < 0.00
Doing a reinit of the equation systems
Advancing to next step
===== Time Step 2 ( 40.00%), time = 0.00 =====
Solving Solid
Assembling the System
Nonlinear Residual: 788389.30
Linear solve starting, tolerance 0.00
Linear solve finished, step 9, residual 0.00
Trying full Newton step
Current Residual: 71.64
Nonlinear step: |du|/|u| = 0.06, |du| = 1.05
Assembling the System
Nonlinear Residual: 71.64
Linear solve starting, tolerance 0.00
Linear solve finished, step 9, residual 0.00
Trying full Newton step
Current Residual: 3.92
Nonlinear solver current_residual 3.92 > 0.00
Nonlinear solver relative residual 0.00 > 0.00
Nonlinear solver converged, step 1, relative step size 0.00 < 0.00
Doing a reinit of the equation systems
Advancing to next step
===== Time Step 3 ( 60.00%), time = 0.00 =====
Solving Solid
Assembling the System
Nonlinear Residual: 844777.18
Linear solve starting, tolerance 0.00
Linear solve finished, step 9, residual 0.00
Trying full Newton step
Current Residual: 120.14
Nonlinear step: |du|/|u| = 0.06, |du| = 1.06
Assembling the System
Nonlinear Residual: 120.14
Linear solve starting, tolerance 0.00
Linear solve finished, step 9, residual 0.00
Trying full Newton step
Current Residual: 8.73
Nonlinear solver current_residual 8.73 > 0.00
Nonlinear solver relative residual 0.00 > 0.00
Nonlinear solver converged, step 1, relative step size 0.00 < 0.00
Doing a reinit of the equation systems
Advancing to next step
===== Time Step 4 ( 80.00%), time = 0.00 =====
Solving Solid
Assembling the System
Nonlinear Residual: 909711.00
Linear solve starting, tolerance 0.00
Linear solve finished, step 9, residual 0.00
Trying full Newton step
Current Residual: 205.60
Nonlinear step: |du|/|u| = 0.07, |du| = 1.07
Assembling the System
Nonlinear Residual: 205.60
Linear solve starting, tolerance 0.00
Linear solve finished, step 9, residual 0.00
Trying full Newton step
Current Residual: 20.81
Nonlinear step: |du|/|u| = 0.00, |du| = 0.02
Assembling the System
Nonlinear Residual: 20.81
Linear solve starting, tolerance 0.00
Linear solve finished, step 10, residual 0.00
Trying full Newton step
Current Residual: 0.14
Nonlinear solver current_residual 0.14 > 0.00
Nonlinear solver relative residual 0.00 > 0.00
Nonlinear solver converged, step 2, relative step size 0.00 < 0.00
Doing a reinit of the equation systems
Advancing to next step
Finished calculations
************************************************************************************************************************
*** WIDEN YOUR WINDOW TO 120 CHARACTERS. Use 'enscript -r -fCourier9' to print this document ***
************************************************************************************************************************
---------------------------------------------- PETSc Performance Summary: ----------------------------------------------
/workspace/libmesh/examples/fem_system/fem_system_ex2/.libs/lt-example-devel on a intel-12. named hbar.ices.utexas.edu with 2 processors, by benkirk Tue Feb 5 13:41:13 2013
Using Petsc Release Version 3.3.0, Patch 2, Fri Jul 13 15:42:00 CDT 2012
Max Max/Min Avg Total
Time (sec): 3.578e+00 1.00000 3.578e+00
Objects: 6.130e+02 1.00000 6.130e+02
Flops: 2.320e+07 3.04748 1.541e+07 3.081e+07
Flops/sec: 6.483e+06 3.04748 4.305e+06 8.611e+06
MPI Messages: 4.100e+02 1.00000 4.100e+02 8.200e+02
MPI Message Lengths: 6.265e+05 1.01953 1.513e+03 1.241e+06
MPI Reductions: 1.665e+03 1.00000
Flop counting convention: 1 flop = 1 real number operation of type (multiply/divide/add/subtract)
e.g., VecAXPY() for real vectors of length N --> 2N flops
and VecAXPY() for complex vectors of length N --> 8N flops
Summary of Stages: ----- Time ------ ----- Flops ----- --- Messages --- -- Message Lengths -- -- Reductions --
Avg %Total Avg %Total counts %Total Avg %Total counts %Total
0: Main Stage: 3.5783e+00 100.0% 3.0812e+07 100.0% 8.200e+02 100.0% 1.513e+03 100.0% 1.664e+03 99.9%
------------------------------------------------------------------------------------------------------------------------
See the 'Profiling' chapter of the users' manual for details on interpreting output.
Phase summary info:
Count: number of times phase was executed
Time and Flops: Max - maximum over all processors
Ratio - ratio of maximum to minimum over all processors
Mess: number of messages sent
Avg. len: average message length
Reduct: number of global reductions
Global: entire computation
Stage: stages of a computation. Set stages with PetscLogStagePush() and PetscLogStagePop().
%T - percent time in this phase %f - percent flops in this phase
%M - percent messages in this phase %L - percent message lengths in this phase
%R - percent reductions in this phase
Total Mflop/s: 10e-6 * (sum of flops over all processors)/(max time over all processors)
------------------------------------------------------------------------------------------------------------------------
Event Count Time (sec) Flops --- Global --- --- Stage --- Total
Max Ratio Max Ratio Max Ratio Mess Avg len Reduct %T %f %M %L %R %T %f %M %L %R Mflop/s
------------------------------------------------------------------------------------------------------------------------
--- Event Stage 0: Main Stage
KSPGMRESOrthog 91 1.0 1.5566e-03 3.4 4.14e+05 1.5 0.0e+00 0.0e+00 9.1e+01 0 2 0 0 5 0 2 0 0 5 443
KSPSetUp 20 1.0 1.3089e-04 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 0
KSPSolve 10 1.0 4.5362e-02 1.0 2.32e+07 3.1 2.0e+02 6.0e+02 2.8e+02 1100 25 10 17 1100 25 10 17 678
PCSetUp 20 1.0 4.0097e-02 3.3 1.48e+07 4.2 0.0e+00 0.0e+00 8.0e+01 1 60 0 0 5 1 60 0 0 5 458
PCSetUpOnBlocks 10 1.0 3.9397e-02 3.5 1.48e+07 4.2 0.0e+00 0.0e+00 7.0e+01 1 60 0 0 4 1 60 0 0 4 467
PCApply 111 1.0 2.5234e-03 1.6 5.37e+06 2.5 0.0e+00 0.0e+00 0.0e+00 0 24 0 0 0 0 24 0 0 0 2988
VecMDot 91 1.0 1.4229e-03 4.8 2.07e+05 1.5 0.0e+00 0.0e+00 9.1e+01 0 1 0 0 5 0 1 0 0 5 242
VecNorm 151 1.0 7.2289e-04 1.8 6.80e+04 1.5 0.0e+00 0.0e+00 1.5e+02 0 0 0 0 9 0 0 0 0 9 157
VecScale 101 1.0 5.7220e-05 1.1 2.27e+04 1.5 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 662
VecCopy 106 1.0 5.5790e-05 1.2 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 0
VecSet 261 1.0 7.8440e-05 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 0
VecAXPY 30 1.0 5.7459e-05 1.9 1.35e+04 1.5 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 392
VecMAXPY 101 1.0 9.2268e-05 1.2 2.48e+05 1.5 0.0e+00 0.0e+00 0.0e+00 0 1 0 0 0 0 1 0 0 0 4479
VecAssemblyBegin 269 1.0 1.1523e-02 1.0 0.00e+00 0.0 6.2e+01 9.2e+02 7.5e+02 0 0 8 5 45 0 0 8 5 45 0
VecAssemblyEnd 269 1.0 1.5426e-04 1.2 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 0
VecScatterBegin 267 1.0 3.6383e-04 1.1 0.00e+00 0.0 4.6e+02 7.5e+02 0.0e+00 0 0 57 28 0 0 0 57 28 0 0
VecScatterEnd 267 1.0 2.8325e-02104.1 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 0
VecNormalize 101 1.0 4.6802e-04 1.4 6.82e+04 1.5 0.0e+00 0.0e+00 1.0e+02 0 0 0 0 6 0 0 0 0 6 243
MatMult 101 1.0 2.8925e-0223.8 2.44e+06 1.6 2.0e+02 6.0e+02 0.0e+00 0 13 25 10 0 0 13 25 10 0 137
MatSolve 111 1.0 1.8594e-03 2.1 5.37e+06 2.5 0.0e+00 0.0e+00 0.0e+00 0 24 0 0 0 0 24 0 0 0 4055
MatLUFactorNum 10 1.0 7.1380e-03 3.3 1.48e+07 4.2 0.0e+00 0.0e+00 0.0e+00 0 60 0 0 0 0 60 0 0 0 2576
MatILUFactorSym 10 1.0 3.1180e-02 3.9 0.00e+00 0.0 0.0e+00 0.0e+00 3.0e+01 1 0 0 0 2 1 0 0 0 2 0
MatAssemblyBegin 20 1.0 2.0442e-03 4.9 0.00e+00 0.0 3.0e+01 2.5e+04 4.0e+01 0 0 4 59 2 0 0 4 59 2 0
MatAssemblyEnd 20 1.0 3.0339e-03 2.0 0.00e+00 0.0 2.0e+01 1.5e+02 4.0e+01 0 0 2 0 2 0 0 2 0 2 0
MatGetRowIJ 10 1.0 2.1458e-06 1.1 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 0
MatGetOrdering 10 1.0 3.9315e-04 1.1 0.00e+00 0.0 0.0e+00 0.0e+00 4.0e+01 0 0 0 0 2 0 0 0 0 2 0
MatZeroEntries 22 1.0 1.8334e-04 2.1 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 0
------------------------------------------------------------------------------------------------------------------------
Memory usage is given in bytes:
Object Type Creations Destructions Memory Descendants' Mem.
Reports information only for process 0.
--- Event Stage 0: Main Stage
Krylov Solver 16 16 104192 0
Preconditioner 16 16 13864 0
Vector 283 283 857384 0
Vector Scatter 71 71 73556 0
Index Set 162 162 150500 0
IS L to G Mapping 36 36 20304 0
Matrix 28 28 4014552 0
Viewer 1 0 0 0
========================================================================================================================
Average time to get PetscTime(): 9.53674e-08
Average time for MPI_Barrier(): 1.23978e-06
Average time for zero size MPI_Send(): 1.20401e-05
#PETSc Option Table entries:
-ksp_right_pc
-log_summary
-pc_type bjacobi
-sub_pc_factor_levels 4
-sub_pc_factor_zeropivot 0
-sub_pc_type ilu
#End of PETSc Option Table entries
Compiled without FORTRAN kernels
Compiled with full precision matrices (default)
sizeof(short) 2 sizeof(int) 4 sizeof(long) 8 sizeof(void*) 8 sizeof(PetscScalar) 8 sizeof(PetscInt) 4
Configure run at: Thu Nov 8 11:21:02 2012
Configure options: --with-debugging=false --COPTFLAGS=-O3 --CXXOPTFLAGS=-O3 --FOPTFLAGS=-O3 --with-clanguage=C++ --with-shared-libraries=1 --with-mpi-dir=/opt/apps/ossw/libraries/mpich2/mpich2-1.4.1p1/sl6/intel-12.1 --with-mumps=true --download-mumps=1 --with-metis=true --download-metis=1 --with-parmetis=true --download-parmetis=1 --with-superlu=true --download-superlu=1 --with-superludir=true --download-superlu_dist=1 --with-blacs=true --download-blacs=1 --with-scalapack=true --download-scalapack=1 --with-hypre=true --download-hypre=1 --with-blas-lib="[/opt/apps/sysnet/intel/12.1/mkl/10.3.12.361/lib/intel64/libmkl_intel_lp64.so,/opt/apps/sysnet/intel/12.1/mkl/10.3.12.361/lib/intel64/libmkl_sequential.so,/opt/apps/sysnet/intel/12.1/mkl/10.3.12.361/lib/intel64/libmkl_core.so]" --with-lapack-lib="[/opt/apps/sysnet/intel/12.1/mkl/10.3.12.361/lib/intel64/libmkl_lapack95_lp64.a]"
-----------------------------------------
Libraries compiled on Thu Nov 8 11:21:02 2012 on daedalus.ices.utexas.edu
Machine characteristics: Linux-2.6.32-279.1.1.el6.x86_64-x86_64-with-redhat-6.3-Carbon
Using PETSc directory: /opt/apps/ossw/libraries/petsc/petsc-3.3-p2
Using PETSc arch: intel-12.1-mkl-intel-10.3.12.361-mpich2-1.4.1p1-cxx-opt
-----------------------------------------
Using C compiler: /opt/apps/ossw/libraries/mpich2/mpich2-1.4.1p1/sl6/intel-12.1/bin/mpicxx -wd1572 -O3 -fPIC ${COPTFLAGS} ${CFLAGS}
Using Fortran compiler: /opt/apps/ossw/libraries/mpich2/mpich2-1.4.1p1/sl6/intel-12.1/bin/mpif90 -fPIC -O3 ${FOPTFLAGS} ${FFLAGS}
-----------------------------------------
Using include paths: -I/opt/apps/ossw/libraries/petsc/petsc-3.3-p2/intel-12.1-mkl-intel-10.3.12.361-mpich2-1.4.1p1-cxx-opt/include -I/opt/apps/ossw/libraries/petsc/petsc-3.3-p2/include -I/opt/apps/ossw/libraries/petsc/petsc-3.3-p2/include -I/opt/apps/ossw/libraries/petsc/petsc-3.3-p2/intel-12.1-mkl-intel-10.3.12.361-mpich2-1.4.1p1-cxx-opt/include -I/opt/apps/ossw/libraries/mpich2/mpich2-1.4.1p1/sl6/intel-12.1/include
-----------------------------------------
Using C linker: /opt/apps/ossw/libraries/mpich2/mpich2-1.4.1p1/sl6/intel-12.1/bin/mpicxx
Using Fortran linker: /opt/apps/ossw/libraries/mpich2/mpich2-1.4.1p1/sl6/intel-12.1/bin/mpif90
Using libraries: -Wl,-rpath,/opt/apps/ossw/libraries/petsc/petsc-3.3-p2/intel-12.1-mkl-intel-10.3.12.361-mpich2-1.4.1p1-cxx-opt/lib -L/opt/apps/ossw/libraries/petsc/petsc-3.3-p2/intel-12.1-mkl-intel-10.3.12.361-mpich2-1.4.1p1-cxx-opt/lib -lpetsc -lX11 -Wl,-rpath,/opt/apps/ossw/libraries/petsc/petsc-3.3-p2/intel-12.1-mkl-intel-10.3.12.361-mpich2-1.4.1p1-cxx-opt/lib -L/opt/apps/ossw/libraries/petsc/petsc-3.3-p2/intel-12.1-mkl-intel-10.3.12.361-mpich2-1.4.1p1-cxx-opt/lib -lcmumps -ldmumps -lsmumps -lzmumps -lmumps_common -lpord -lHYPRE -lpthread -lsuperlu_dist_3.0 -lparmetis -lmetis -lscalapack -lblacs -lsuperlu_4.3 -Wl,-rpath,/opt/apps/sysnet/intel/12.1/mkl/10.3.12.361/lib/intel64 -L/opt/apps/sysnet/intel/12.1/mkl/10.3.12.361/lib/intel64 -lmkl_lapack95_lp64 -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -Wl,-rpath,/opt/apps/ossw/libraries/mpich2/mpich2-1.4.1p1/sl6/intel-12.1/lib -L/opt/apps/ossw/libraries/mpich2/mpich2-1.4.1p1/sl6/intel-12.1/lib -Wl,-rpath,/opt/apps/sysnet/intel/12.1/composer_xe_2011_sp1.7.256/compiler/lib/intel64 -L/opt/apps/sysnet/intel/12.1/composer_xe_2011_sp1.7.256/compiler/lib/intel64 -Wl,-rpath,/usr/lib/gcc/x86_64-redhat-linux/4.4.6 -L/usr/lib/gcc/x86_64-redhat-linux/4.4.6 -lmpichf90 -lifport -lifcore -lm -lm -lmpichcxx -ldl -lmpich -lopa -lmpl -lrt -lpthread -limf -lsvml -lipgo -ldecimal -lcilkrts -lstdc++ -lgcc_s -lirc -lirc_s -ldl
-----------------------------------------
----------------------------------------------------------------------------------------------------------------------
| Processor id: 0 |
| Num Processors: 2 |
| Time: Tue Feb 5 13:41:13 2013 |
| OS: Linux |
| HostName: hbar.ices.utexas.edu |
| OS Release: 2.6.32-279.1.1.el6.x86_64 |
| OS Version: #1 SMP Tue Jul 10 11:24:23 CDT 2012 |
| Machine: x86_64 |
| Username: benkirk |
| Configuration: ./configure '--enable-everything' |
| '--prefix=/workspace/libmesh/install' |
| 'CXX=mpicxx' |
| 'CC=mpicc' |
| 'F77=mpif77' |
| 'FC=mpif90' |
| 'PETSC_DIR=/opt/apps/ossw/libraries/petsc/petsc-3.3-p2' |
| 'PETSC_ARCH=intel-12.1-mkl-intel-10.3.12.361-mpich2-1.4.1p1-cxx-opt' |
| 'SLEPC_DIR=/opt/apps/ossw/libraries/slepc/slepc-3.3-p2-petsc-3.3-p2-cxx-opt' |
| 'TRILINOS_DIR=/opt/apps/ossw/libraries/trilinos/trilinos-10.12.2/sl6/intel-12.1/mpich2-1.4.1p1/mkl-intel-10.3.12.361'|
| 'VTK_DIR=/opt/apps/ossw/libraries/vtk/vtk-5.10.0/sl6/intel-12.1' |
----------------------------------------------------------------------------------------------------------------------
----------------------------------------------------------------------------------------------------------------
| libMesh Performance: Alive time=3.58831, Active time=3.54079 |
----------------------------------------------------------------------------------------------------------------
| Event nCalls Total Time Avg Time Total Time Avg Time % of Active Time |
| w/o Sub w/o Sub With Sub With Sub w/o S With S |
|----------------------------------------------------------------------------------------------------------------|
| |
| |
| DofMap |
| add_neighbors_to_send_list() 12 0.0497 0.004138 0.1113 0.009273 1.40 3.14 |
| build_sparsity() 6 0.0451 0.007512 0.0862 0.014365 1.27 2.43 |
| create_dof_constraints() 12 0.0018 0.000153 0.0018 0.000153 0.05 0.05 |
| distribute_dofs() 12 0.0170 0.001416 0.0436 0.003634 0.48 1.23 |
| dof_indices() 27584 1.7371 0.000063 1.7371 0.000063 49.06 49.06 |
| old_dof_indices() 2560 0.2120 0.000083 0.2120 0.000083 5.99 5.99 |
| prepare_send_list() 12 0.0011 0.000088 0.0011 0.000088 0.03 0.03 |
| reinit() 12 0.0241 0.002004 0.0241 0.002004 0.68 0.68 |
| |
| EquationSystems |
| build_solution_vector() 5 0.0042 0.000849 0.0529 0.010585 0.12 1.49 |
| |
| FE |
| compute_shape_functions() 5760 0.3902 0.000068 0.3902 0.000068 11.02 11.02 |
| init_shape_functions() 2040 0.0470 0.000023 0.0470 0.000023 1.33 1.33 |
| |
| FEMSystem |
| assembly() 10 0.5345 0.053451 1.4211 0.142108 15.10 40.13 |
| assembly(get_residual) 10 0.1158 0.011577 1.0105 0.101051 3.27 28.54 |
| |
| FEMap |
| compute_affine_map() 290 0.0044 0.000015 0.0044 0.000015 0.12 0.12 |
| compute_face_map() 1920 0.0227 0.000012 0.0227 0.000012 0.64 0.64 |
| compute_map() 5470 0.1138 0.000021 0.1138 0.000021 3.21 3.21 |
| init_face_shape_functions() 40 0.0010 0.000025 0.0010 0.000025 0.03 0.03 |
| init_reference_to_physical_map() 2040 0.0667 0.000033 0.0667 0.000033 1.88 1.88 |
| |
| LocationMap |
| init() 5 0.0006 0.000129 0.0006 0.000129 0.02 0.02 |
| |
| Mesh |
| contract() 5 0.0001 0.000021 0.0003 0.000054 0.00 0.01 |
| find_neighbors() 1 0.0017 0.001743 0.0018 0.001810 0.05 0.05 |
| renumber_nodes_and_elem() 7 0.0003 0.000043 0.0003 0.000043 0.01 0.01 |
| |
| MeshCommunication |
| compute_hilbert_indices() 2 0.0007 0.000329 0.0007 0.000329 0.02 0.02 |
| find_global_indices() 2 0.0005 0.000243 0.0021 0.001041 0.01 0.06 |
| parallel_sort() 2 0.0006 0.000285 0.0007 0.000326 0.02 0.02 |
| |
| MeshOutput |
| write_equation_systems() 5 0.0104 0.002070 0.0635 0.012709 0.29 1.79 |
| |
| MeshRefinement |
| _coarsen_elements() 5 0.0001 0.000018 0.0001 0.000022 0.00 0.00 |
| _refine_elements() 5 0.0002 0.000036 0.0002 0.000045 0.01 0.01 |
| make_coarsening_compatible() 10 0.0003 0.000025 0.0003 0.000025 0.01 0.01 |
| make_refinement_compatible() 10 0.0000 0.000002 0.0000 0.000004 0.00 0.00 |
| |
| MeshTools::Generation |
| build_cube() 1 0.0007 0.000673 0.0007 0.000673 0.02 0.02 |
| |
| MetisPartitioner |
| partition() 1 0.0014 0.001419 0.0021 0.002147 0.04 0.06 |
| |
| NewtonSolver |
| solve() 5 0.0283 0.005668 2.7814 0.556278 0.80 78.55 |
| |
| Parallel |
| allgather() 40 0.0002 0.000006 0.0003 0.000006 0.01 0.01 |
| max(bool) 21 0.0001 0.000005 0.0001 0.000005 0.00 0.00 |
| max(scalar) 1165 0.0030 0.000003 0.0030 0.000003 0.08 0.08 |
| max(vector) 283 0.0018 0.000006 0.0040 0.000014 0.05 0.11 |
| min(bool) 1453 0.0036 0.000002 0.0036 0.000002 0.10 0.10 |
| min(scalar) 1152 0.0199 0.000017 0.0199 0.000017 0.56 0.56 |
| min(vector) 283 0.0020 0.000007 0.0043 0.000015 0.06 0.12 |
| probe() 146 0.0014 0.000010 0.0014 0.000010 0.04 0.04 |
| receive() 146 0.0008 0.000005 0.0022 0.000015 0.02 0.06 |
| send() 146 0.0004 0.000003 0.0004 0.000003 0.01 0.01 |
| send_receive() 150 0.0012 0.000008 0.0040 0.000027 0.03 0.11 |
| sum() 66 0.0004 0.000006 0.0006 0.000010 0.01 0.02 |
| |
| Parallel::Request |
| wait() 146 0.0002 0.000002 0.0002 0.000002 0.01 0.01 |
| |
| Partitioner |
| set_node_processor_ids() 1 0.0003 0.000330 0.0004 0.000410 0.01 0.01 |
| set_parent_processor_ids() 1 0.0001 0.000092 0.0001 0.000092 0.00 0.00 |
| |
| PetscLinearSolver |
| solve() 10 0.0497 0.004966 0.0497 0.004966 1.40 1.40 |
| |
| ProjectVector |
| operator() 20 0.0136 0.000678 0.2300 0.011498 0.38 6.49 |
| |
| System |
| project_vector() 20 0.0084 0.000419 0.3446 0.017231 0.24 9.73 |
----------------------------------------------------------------------------------------------------------------
| Totals: 53110 3.5408 100.00 |
----------------------------------------------------------------------------------------------------------------
***************************************************************
* Done Running Example fem_system_ex2:
* mpirun -np 2 example-devel -pc_type bjacobi -sub_pc_type ilu -sub_pc_factor_levels 4 -sub_pc_factor_zeropivot 0 -ksp_right_pc -log_summary
***************************************************************