The source file transient_ex2.C with comments:
#include <iostream>
#include <fstream>
#include <algorithm>
#include <stdio.h>
#include <math.h>
Basic include file needed for the mesh functionality.
#include "libmesh/libmesh.h"
#include "libmesh/serial_mesh.h"
#include "libmesh/gmv_io.h"
#include "libmesh/vtk_io.h"
#include "libmesh/newmark_system.h"
#include "libmesh/equation_systems.h"
Define the Finite Element object.
#include "libmesh/fe.h"
Define Gauss quadrature rules.
#include "libmesh/quadrature_gauss.h"
Define useful datatypes for finite element
#include "libmesh/dense_matrix.h"
#include "libmesh/dense_vector.h"
Define matrix and vector data types for the global
equation system. These are base classes,
from which specific implementations, like
the PETSc or LASPACK implementations, are derived.
#include "libmesh/sparse_matrix.h"
#include "libmesh/numeric_vector.h"
Define the DofMap, which handles degree of freedom
indexing.
#include "libmesh/dof_map.h"
The definition of a vertex associated with a Mesh.
#include "libmesh/node.h"
The definition of a geometric element
#include "libmesh/elem.h"
Defines the MeshData class, which allows you to store
data about the mesh when reading in files, etc.
#include "libmesh/mesh_data.h"
Bring in everything from the libMesh namespace
using namespace libMesh;
Function prototype. This is the function that will assemble
the linear system for our problem, governed by the linear
wave equation.
void assemble_wave(EquationSystems& es,
const std::string& system_name);
Function Prototype. This function will be used to apply the
initial conditions.
void apply_initial(EquationSystems& es,
const std::string& system_name);
Function Prototype. This function imposes
Dirichlet Boundary conditions via the penalty
method after the system is assembled.
void fill_dirichlet_bc(EquationSystems& es,
const std::string& system_name);
The main program
int main (int argc, char** argv)
{
Initialize libraries, like in example 2.
LibMeshInit init (argc, argv);
Check for proper usage.
if (argc < 2)
{
if (libMesh::processor_id() == 0)
std::cerr << "Usage: " << argv[0] << " [meshfile]"
<< std::endl;
libmesh_error();
}
Tell the user what we are doing.
else
{
std::cout << "Running " << argv[0];
for (int i=1; i<argc; i++)
std::cout << " " << argv[i];
std::cout << std::endl << std::endl;
}
LasPack solvers don't work so well for this example
(not sure why), and Trilinos matrices don't work at all.
Print a warning to the user if PETSc is not in use.
if (libMesh::default_solver_package() == LASPACK_SOLVERS)
{
std::cout << "WARNING! It appears you are using the\n"
<< "LasPack solvers. This example may not converge\n"
<< "using LasPack, but should work OK with PETSc.\n"
<< "http://www.mcs.anl.gov/petsc/\n"
<< std::endl;
}
else if (libMesh::default_solver_package() == TRILINOS_SOLVERS)
{
std::cout << "WARNING! It appears you are using the\n"
<< "Trilinos solvers. The current libMesh Epetra\n"
<< "interface does not allow sparse matrix addition,\n"
<< "as is needed in this problem. We recommend\n"
<< "using PETSc: http://www.mcs.anl.gov/petsc/\n"
<< std::endl;
return 0;
}
Get the name of the mesh file
from the command line.
std::string mesh_file = argv[1];
std::cout << "Mesh file is: " << mesh_file << std::endl;
Skip this 3D example if libMesh was compiled as 1D or 2D-only.
libmesh_example_assert(3 <= LIBMESH_DIM, "3D support");
Create a mesh.
This example directly references all mesh nodes and is
incompatible with ParallelMesh use.
SerialMesh mesh;
MeshData mesh_data(mesh);
Read the meshfile specified in the command line or
use the internal mesh generator to create a uniform
grid on an elongated cube.
mesh.read(mesh_file, &mesh_data);
mesh.build_cube (10, 10, 40,
-1., 1.,
-1., 1.,
0., 4.,
HEX8);
Print information about the mesh to the screen.
Print information about the mesh to the screen.
mesh.print_info();
The node that should be monitored.
const unsigned int result_node = 274;
Time stepping issues
Note that the total current time is stored as a parameter in the \pEquationSystems object.
the time step size
Note that the total current time is stored as a parameter in the \pEquationSystems object.
the time step size
const Real delta_t = .0000625;
The number of time steps.
unsigned int n_time_steps = 300;
Create an equation systems object.
EquationSystems equation_systems (mesh);
Declare the system and its variables.
Create a NewmarkSystem named "Wave"
equation_systems.add_system<NewmarkSystem> ("Wave");
Use a handy reference to this system
NewmarkSystem & t_system = equation_systems.get_system<NewmarkSystem> ("Wave");
Add the variable "p" to "Wave". "p"
will be approximated using first-order approximation.
t_system.add_variable("p", FIRST);
Give the system a pointer to the matrix assembly
function and the initial condition function defined
below.
t_system.attach_assemble_function (assemble_wave);
t_system.attach_init_function (apply_initial);
Set the time step size, and optionally the
Newmark parameters, so that \p NewmarkSystem can
compute integration constants. Here we simply use
pass only the time step and use default values
for \p alpha=.25 and \p delta=.5.
t_system.set_newmark_parameters(delta_t);
Set the speed of sound and fluid density
as \p EquationSystems parameter,
so that \p assemble_wave() can access it.
equation_systems.parameters.set<Real>("speed") = 1000.;
equation_systems.parameters.set<Real>("fluid density") = 1000.;
Start time integration from t=0
t_system.time = 0.;
Initialize the data structures for the equation system.
equation_systems.init();
Prints information about the system to the screen.
equation_systems.print_info();
A file to store the results at certain nodes.
std::ofstream res_out("pressure_node.res");
get the dof_numbers for the nodes that
should be monitored.
const unsigned int res_node_no = result_node;
const Node& res_node = mesh.node(res_node_no-1);
unsigned int dof_no = res_node.dof_number(0,0,0);
Assemble the time independent system matrices and rhs.
This function will also compute the effective system matrix
K~=K+a_0*M+a_1*C and apply user specified initial
conditions.
t_system.assemble();
Now solve for each time step.
For convenience, use a local buffer of the
current time. But once this time is updated,
also update the \p EquationSystems parameter
Start with t_time = 0 and write a short header
to the nodal result file
res_out << "# pressure at node " << res_node_no << "\n"
<< "# time\tpressure\n"
<< t_system.time << "\t" << 0 << std::endl;
for (unsigned int time_step=0; time_step<n_time_steps; time_step++)
{
Update the time. Both here and in the
\p EquationSystems object
t_system.time += delta_t;
Update the rhs.
t_system.update_rhs();
Impose essential boundary conditions.
Not that since the matrix is only assembled once,
the penalty parameter should be added to the matrix
only in the first time step. The applied
boundary conditions may be time-dependent and hence
the rhs vector is considered in each time step.
if (time_step == 0)
{
The local function \p fill_dirichlet_bc()
may also set Dirichlet boundary conditions for the
matrix. When you set the flag as shown below,
the flag will return true. If you want it to return
false, simply do not set it.
equation_systems.parameters.set<bool>("Newmark set BC for Matrix") = true;
fill_dirichlet_bc(equation_systems, "Wave");
unset the flag, so that it returns false
equation_systems.parameters.set<bool>("Newmark set BC for Matrix") = false;
}
else
fill_dirichlet_bc(equation_systems, "Wave");
Solve the system "Wave".
t_system.solve();
After solving the system, write the solution
to a GMV-formatted plot file.
Do only for a few time steps.
if (time_step == 30 || time_step == 60 ||
time_step == 90 || time_step == 120 )
{
char buf[14];
if (!libMesh::on_command_line("--vtk")){
sprintf (buf, "out.%03d.gmv", time_step);
GMVIO(mesh).write_equation_systems (buf,equation_systems);
}else{
#ifdef LIBMESH_HAVE_VTK
VTK viewers are generally not happy with two dots in a filename
sprintf (buf, "out_%03d.exd", time_step);
VTKIO(mesh).write_equation_systems (buf,equation_systems);
#endif // #ifdef LIBMESH_HAVE_VTK
}
}
Update the p, v and a.
t_system.update_u_v_a();
dof_no may not be local in parallel runs, so we may need a
global displacement vector
NumericVector<Number> &displacement
= t_system.get_vector("displacement");
std::vector<Number> global_displacement(displacement.size());
displacement.localize(global_displacement);
Write nodal results to file. The results can then
be viewed with e.g. gnuplot (run gnuplot and type
'plot "pressure_node.res" with lines' in the command line)
res_out << t_system.time << "\t"
<< global_displacement[dof_no]
<< std::endl;
}
All done.
return 0;
}
This function assembles the system matrix and right-hand-side
for our wave equation.
void assemble_wave(EquationSystems& es,
const std::string& system_name)
{
It is a good idea to make sure we are assembling
the proper system.
libmesh_assert_equal_to (system_name, "Wave");
Get a constant reference to the mesh object.
const MeshBase& mesh = es.get_mesh();
The dimension that we are running.
const unsigned int dim = mesh.mesh_dimension();
Copy the speed of sound and fluid density
to a local variable.
const Real speed = es.parameters.get<Real>("speed");
const Real rho = es.parameters.get<Real>("fluid density");
Get a reference to our system, as before.
NewmarkSystem & t_system = es.get_system<NewmarkSystem> (system_name);
Get a constant reference to the Finite Element type
for the first (and only) variable in the system.
FEType fe_type = t_system.get_dof_map().variable_type(0);
In here, we will add the element matrices to the
@e additional matrices "stiffness_mass" and "damping"
and the additional vector "force", not to the members
"matrix" and "rhs". Therefore, get writable
references to them.
SparseMatrix<Number>& stiffness = t_system.get_matrix("stiffness");
SparseMatrix<Number>& damping = t_system.get_matrix("damping");
SparseMatrix<Number>& mass = t_system.get_matrix("mass");
NumericVector<Number>& force = t_system.get_vector("force");
Some solver packages (PETSc) are especially picky about
allocating sparsity structure and truly assigning values
to this structure. Namely, matrix additions, as performed
later, exhibit acceptable performance only for identical
sparsity structures. Therefore, explicitly zero the
values in the collective matrix, so that matrix additions
encounter identical sparsity structures.
SparseMatrix<Number>& matrix = *t_system.matrix;
DenseMatrix<Number> zero_matrix;
Build a Finite Element object of the specified type. Since the
\p FEBase::build() member dynamically creates memory we will
store the object as an \p AutoPtr. This can be thought
of as a pointer that will clean up after itself.
AutoPtr<FEBase> fe (FEBase::build(dim, fe_type));
A 2nd order Gauss quadrature rule for numerical integration.
QGauss qrule (dim, SECOND);
Tell the finite element object to use our quadrature rule.
fe->attach_quadrature_rule (&qrule);
The element Jacobian * quadrature weight at each integration point.
const std::vector<Real>& JxW = fe->get_JxW();
The element shape functions evaluated at the quadrature points.
const std::vector<std::vector<Real> >& phi = fe->get_phi();
The element shape function gradients evaluated at the quadrature
points.
const std::vector<std::vector<RealGradient> >& dphi = fe->get_dphi();
A reference to the \p DofMap object for this system. The \p DofMap
object handles the index translation from node and element numbers
to degree of freedom numbers.
const DofMap& dof_map = t_system.get_dof_map();
The element mass, damping and stiffness matrices
and the element contribution to the rhs.
DenseMatrix<Number> Ke, Ce, Me;
DenseVector<Number> Fe;
This vector will hold the degree of freedom indices for
the element. These define where in the global system
the element degrees of freedom get mapped.
std::vector<dof_id_type> dof_indices;
Now we will loop over all the elements in the mesh.
We will compute the element matrix and right-hand-side
contribution.
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)
{
Store a pointer to the element we are currently
working on. This allows for nicer syntax later.
const Elem* elem = *el;
Get the degree of freedom indices for the
current element. These define where in the global
matrix and right-hand-side this element will
contribute to.
dof_map.dof_indices (elem, dof_indices);
Compute the element-specific data for the current
element. This involves computing the location of the
quadrature points (q_point) and the shape functions
(phi, dphi) for the current element.
fe->reinit (elem);
Zero the element matrices and rhs before
summing them. We use the resize member here because
the number of degrees of freedom might have changed from
the last element. Note that this will be the case if the
element type is different (i.e. the last element was HEX8
and now have a PRISM6).
{
const unsigned int n_dof_indices = dof_indices.size();
Ke.resize (n_dof_indices, n_dof_indices);
Ce.resize (n_dof_indices, n_dof_indices);
Me.resize (n_dof_indices, n_dof_indices);
zero_matrix.resize (n_dof_indices, n_dof_indices);
Fe.resize (n_dof_indices);
}
Now loop over the quadrature points. This handles
the numeric integration.
for (unsigned int qp=0; qp<qrule.n_points(); qp++)
{
Now we will build the element matrix. This involves
a double loop to integrate the test funcions (i) against
the trial functions (j).
for (unsigned int i=0; i<phi.size(); i++)
for (unsigned int j=0; j<phi.size(); j++)
{
Ke(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
Me(i,j) += JxW[qp]*phi[i][qp]*phi[j][qp]
*1./(speed*speed);
} // end of the matrix summation loop
} // end of quadrature point loop
Now compute the contribution to the element matrix and the
right-hand-side vector if the current element lies on the
boundary.
{
In this example no natural boundary conditions will
be considered. The code is left here so it can easily
be extended.
don't do this for any side
don't do this for any side
for (unsigned int side=0; side<elem->n_sides(); side++)
if (!true)
if (elem->neighbor(side) == NULL)
{
Declare a special finite element object for
boundary integration.
AutoPtr<FEBase> fe_face (FEBase::build(dim, fe_type));
Boundary integration requires one quadraure rule,
with dimensionality one less than the dimensionality
of the element.
QGauss qface(dim-1, SECOND);
Tell the finte element object to use our
quadrature rule.
fe_face->attach_quadrature_rule (&qface);
The value of the shape functions at the quadrature
points.
const std::vector<std::vector<Real> >& phi_face = fe_face->get_phi();
The Jacobian * Quadrature Weight at the quadrature
points on the face.
const std::vector<Real>& JxW_face = fe_face->get_JxW();
Compute the shape function values on the element
face.
fe_face->reinit(elem, side);
Here we consider a normal acceleration acc_n=1 applied to
the whole boundary of our mesh.
const Real acc_n_value = 1.0;
Loop over the face quadrature points for integration.
for (unsigned int qp=0; qp<qface.n_points(); qp++)
{
Right-hand-side contribution due to prescribed
normal acceleration.
for (unsigned int i=0; i<phi_face.size(); i++)
{
Fe(i) += acc_n_value*rho
*phi_face[i][qp]*JxW_face[qp];
}
} // end face quadrature point loop
} // end if (elem->neighbor(side) == NULL)
In this example the Dirichlet boundary conditions will be
imposed via panalty method after the
system is assembled.
} // end boundary condition section
If this assembly program were to be used on an adaptive mesh,
we would have to apply any hanging node constraint equations
by uncommenting the following lines:
std::vector dof_indicesC = dof_indices;
std::vector dof_indicesM = dof_indices;
dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices);
dof_map.constrain_element_matrix (Ce, dof_indicesC);
dof_map.constrain_element_matrix (Me, dof_indicesM);
Finally, simply add the contributions to the additional matrices and vector.
Finally, simply add the contributions to the additional matrices and vector.
stiffness.add_matrix (Ke, dof_indices);
damping.add_matrix (Ce, dof_indices);
mass.add_matrix (Me, dof_indices);
force.add_vector (Fe, dof_indices);
For the overall matrix, explicitly zero the entries where
we added values in the other ones, so that we have
identical sparsity footprints.
matrix.add_matrix(zero_matrix, dof_indices);
} // end of element loop
All done!
return;
}
This function applies the initial conditions
void apply_initial(EquationSystems& es,
const std::string& system_name)
{
Get a reference to our system, as before
NewmarkSystem & t_system = es.get_system<NewmarkSystem> (system_name);
Numeric vectors for the pressure, velocity and acceleration
values.
NumericVector<Number>& pres_vec = t_system.get_vector("displacement");
NumericVector<Number>& vel_vec = t_system.get_vector("velocity");
NumericVector<Number>& acc_vec = t_system.get_vector("acceleration");
Assume our fluid to be at rest, which would
also be the default conditions in class NewmarkSystem,
but let us do it explicetly here.
pres_vec.zero();
vel_vec.zero();
acc_vec.zero();
}
This function applies the Dirichlet boundary conditions
void fill_dirichlet_bc(EquationSystems& es,
const std::string& system_name)
{
It is a good idea to make sure we are assembling
the proper system.
libmesh_assert_equal_to (system_name, "Wave");
Get a reference to our system, as before.
NewmarkSystem & t_system = es.get_system<NewmarkSystem> (system_name);
Get writable references to the overall matrix and vector.
SparseMatrix<Number>& matrix = *t_system.matrix;
NumericVector<Number>& rhs = *t_system.rhs;
Get a constant reference to the mesh object.
const MeshBase& mesh = es.get_mesh();
Get \p libMesh's \f$ \pi \f$
const Real pi = libMesh::pi;
Ask the \p EquationSystems flag whether
we should do this also for the matrix
const bool do_for_matrix =
es.parameters.get<bool>("Newmark set BC for Matrix");
Number of nodes in the mesh.
unsigned int n_nodes = mesh.n_nodes();
for (unsigned int n_cnt=0; n_cnt<n_nodes; n_cnt++)
{
Get a reference to the current node.
const Node& curr_node = mesh.node(n_cnt);
Check if Dirichlet BCs should be applied to this node.
Use the \p TOLERANCE from \p mesh_common.h as tolerance.
Here a pressure value is applied if the z-coord.
is equal to 4, which corresponds to one end of the
pipe-mesh in this directory.
const Real z_coo = 4.;
if (fabs(curr_node(2)-z_coo) < TOLERANCE)
{
The global number of the respective degree of freedom.
unsigned int dn = curr_node.dof_number(0,0,0);
The penalty parameter.
const Real penalty = 1.e10;
Here we apply sinusoidal pressure values for 0
Real p_value;
if (t_system.time < .002 )
p_value = sin(2*pi*t_system.time/.002);
else
p_value = .0;
Now add the contributions to the matrix and the rhs.
rhs.add(dn, p_value*penalty);
Add the panalty parameter to the global matrix
if desired.
if (do_for_matrix)
matrix.add(dn, dn, penalty);
}
} // loop n_cnt
}
The source file transient_ex2.C without comments:
#include <iostream>
#include <fstream>
#include <algorithm>
#include <stdio.h>
#include <math.h>
#include "libmesh/libmesh.h"
#include "libmesh/serial_mesh.h"
#include "libmesh/gmv_io.h"
#include "libmesh/vtk_io.h"
#include "libmesh/newmark_system.h"
#include "libmesh/equation_systems.h"
#include "libmesh/fe.h"
#include "libmesh/quadrature_gauss.h"
#include "libmesh/dense_matrix.h"
#include "libmesh/dense_vector.h"
#include "libmesh/sparse_matrix.h"
#include "libmesh/numeric_vector.h"
#include "libmesh/dof_map.h"
#include "libmesh/node.h"
#include "libmesh/elem.h"
#include "libmesh/mesh_data.h"
using namespace libMesh;
void assemble_wave(EquationSystems& es,
const std::string& system_name);
void apply_initial(EquationSystems& es,
const std::string& system_name);
void fill_dirichlet_bc(EquationSystems& es,
const std::string& system_name);
int main (int argc, char** argv)
{
LibMeshInit init (argc, argv);
if (argc < 2)
{
if (libMesh::processor_id() == 0)
std::cerr << "Usage: " << argv[0] << " [meshfile]"
<< std::endl;
libmesh_error();
}
else
{
std::cout << "Running " << argv[0];
for (int i=1; i<argc; i++)
std::cout << " " << argv[i];
std::cout << std::endl << std::endl;
}
if (libMesh::default_solver_package() == LASPACK_SOLVERS)
{
std::cout << "WARNING! It appears you are using the\n"
<< "LasPack solvers. This example may not converge\n"
<< "using LasPack, but should work OK with PETSc.\n"
<< "http://www.mcs.anl.gov/petsc/\n"
<< std::endl;
}
else if (libMesh::default_solver_package() == TRILINOS_SOLVERS)
{
std::cout << "WARNING! It appears you are using the\n"
<< "Trilinos solvers. The current libMesh Epetra\n"
<< "interface does not allow sparse matrix addition,\n"
<< "as is needed in this problem. We recommend\n"
<< "using PETSc: http://www.mcs.anl.gov/petsc/\n"
<< std::endl;
return 0;
}
std::string mesh_file = argv[1];
std::cout << "Mesh file is: " << mesh_file << std::endl;
libmesh_example_assert(3 <= LIBMESH_DIM, "3D support");
SerialMesh mesh;
MeshData mesh_data(mesh);
mesh.read(mesh_file, &mesh_data);
mesh.print_info();
const unsigned int result_node = 274;
const Real delta_t = .0000625;
unsigned int n_time_steps = 300;
EquationSystems equation_systems (mesh);
equation_systems.add_system<NewmarkSystem> ("Wave");
NewmarkSystem & t_system = equation_systems.get_system<NewmarkSystem> ("Wave");
t_system.add_variable("p", FIRST);
t_system.attach_assemble_function (assemble_wave);
t_system.attach_init_function (apply_initial);
t_system.set_newmark_parameters(delta_t);
equation_systems.parameters.set<Real>("speed") = 1000.;
equation_systems.parameters.set<Real>("fluid density") = 1000.;
t_system.time = 0.;
equation_systems.init();
equation_systems.print_info();
std::ofstream res_out("pressure_node.res");
const unsigned int res_node_no = result_node;
const Node& res_node = mesh.node(res_node_no-1);
unsigned int dof_no = res_node.dof_number(0,0,0);
t_system.assemble();
res_out << "# pressure at node " << res_node_no << "\n"
<< "# time\tpressure\n"
<< t_system.time << "\t" << 0 << std::endl;
for (unsigned int time_step=0; time_step<n_time_steps; time_step++)
{
t_system.time += delta_t;
t_system.update_rhs();
if (time_step == 0)
{
equation_systems.parameters.set<bool>("Newmark set BC for Matrix") = true;
fill_dirichlet_bc(equation_systems, "Wave");
equation_systems.parameters.set<bool>("Newmark set BC for Matrix") = false;
}
else
fill_dirichlet_bc(equation_systems, "Wave");
t_system.solve();
if (time_step == 30 || time_step == 60 ||
time_step == 90 || time_step == 120 )
{
char buf[14];
if (!libMesh::on_command_line("--vtk")){
sprintf (buf, "out.%03d.gmv", time_step);
GMVIO(mesh).write_equation_systems (buf,equation_systems);
}else{
#ifdef LIBMESH_HAVE_VTK
sprintf (buf, "out_%03d.exd", time_step);
VTKIO(mesh).write_equation_systems (buf,equation_systems);
#endif // #ifdef LIBMESH_HAVE_VTK
}
}
t_system.update_u_v_a();
NumericVector<Number> &displacement
= t_system.get_vector("displacement");
std::vector<Number> global_displacement(displacement.size());
displacement.localize(global_displacement);
res_out << t_system.time << "\t"
<< global_displacement[dof_no]
<< std::endl;
}
return 0;
}
void assemble_wave(EquationSystems& es,
const std::string& system_name)
{
libmesh_assert_equal_to (system_name, "Wave");
const MeshBase& mesh = es.get_mesh();
const unsigned int dim = mesh.mesh_dimension();
const Real speed = es.parameters.get<Real>("speed");
const Real rho = es.parameters.get<Real>("fluid density");
NewmarkSystem & t_system = es.get_system<NewmarkSystem> (system_name);
FEType fe_type = t_system.get_dof_map().variable_type(0);
SparseMatrix<Number>& stiffness = t_system.get_matrix("stiffness");
SparseMatrix<Number>& damping = t_system.get_matrix("damping");
SparseMatrix<Number>& mass = t_system.get_matrix("mass");
NumericVector<Number>& force = t_system.get_vector("force");
SparseMatrix<Number>& matrix = *t_system.matrix;
DenseMatrix<Number> zero_matrix;
AutoPtr<FEBase> fe (FEBase::build(dim, fe_type));
QGauss qrule (dim, SECOND);
fe->attach_quadrature_rule (&qrule);
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 DofMap& dof_map = t_system.get_dof_map();
DenseMatrix<Number> Ke, Ce, Me;
DenseVector<Number> Fe;
std::vector<dof_id_type> dof_indices;
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);
fe->reinit (elem);
{
const unsigned int n_dof_indices = dof_indices.size();
Ke.resize (n_dof_indices, n_dof_indices);
Ce.resize (n_dof_indices, n_dof_indices);
Me.resize (n_dof_indices, n_dof_indices);
zero_matrix.resize (n_dof_indices, n_dof_indices);
Fe.resize (n_dof_indices);
}
for (unsigned int qp=0; qp<qrule.n_points(); qp++)
{
for (unsigned int i=0; i<phi.size(); i++)
for (unsigned int j=0; j<phi.size(); j++)
{
Ke(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
Me(i,j) += JxW[qp]*phi[i][qp]*phi[j][qp]
*1./(speed*speed);
} // end of the matrix summation loop
} // end of quadrature point loop
{
for (unsigned int side=0; side<elem->n_sides(); side++)
if (!true)
{
AutoPtr<FEBase> fe_face (FEBase::build(dim, fe_type));
QGauss qface(dim-1, SECOND);
fe_face->attach_quadrature_rule (&qface);
const std::vector<std::vector<Real> >& phi_face = fe_face->get_phi();
const std::vector<Real>& JxW_face = fe_face->get_JxW();
fe_face->reinit(elem, side);
const Real acc_n_value = 1.0;
for (unsigned int qp=0; qp<qface.n_points(); qp++)
{
for (unsigned int i=0; i<phi_face.size(); i++)
{
Fe(i) += acc_n_value*rho
*phi_face[i][qp]*JxW_face[qp];
}
} // end face quadrature point loop
} // end if (elem->neighbor(side) == NULL)
} // end boundary condition section
stiffness.add_matrix (Ke, dof_indices);
damping.add_matrix (Ce, dof_indices);
mass.add_matrix (Me, dof_indices);
force.add_vector (Fe, dof_indices);
matrix.add_matrix(zero_matrix, dof_indices);
} // end of element loop
return;
}
void apply_initial(EquationSystems& es,
const std::string& system_name)
{
NewmarkSystem & t_system = es.get_system<NewmarkSystem> (system_name);
NumericVector<Number>& pres_vec = t_system.get_vector("displacement");
NumericVector<Number>& vel_vec = t_system.get_vector("velocity");
NumericVector<Number>& acc_vec = t_system.get_vector("acceleration");
pres_vec.zero();
vel_vec.zero();
acc_vec.zero();
}
void fill_dirichlet_bc(EquationSystems& es,
const std::string& system_name)
{
libmesh_assert_equal_to (system_name, "Wave");
NewmarkSystem & t_system = es.get_system<NewmarkSystem> (system_name);
SparseMatrix<Number>& matrix = *t_system.matrix;
NumericVector<Number>& rhs = *t_system.rhs;
const MeshBase& mesh = es.get_mesh();
const Real pi = libMesh::pi;
const bool do_for_matrix =
es.parameters.get<bool>("Newmark set BC for Matrix");
unsigned int n_nodes = mesh.n_nodes();
for (unsigned int n_cnt=0; n_cnt<n_nodes; n_cnt++)
{
const Node& curr_node = mesh.node(n_cnt);
const Real z_coo = 4.;
if (fabs(curr_node(2)-z_coo) < TOLERANCE)
{
unsigned int dn = curr_node.dof_number(0,0,0);
const Real penalty = 1.e10;
Real p_value;
if (t_system.time < .002 )
p_value = sin(2*pi*t_system.time/.002);
else
p_value = .0;
rhs.add(dn, p_value*penalty);
if (do_for_matrix)
matrix.add(dn, dn, penalty);
}
} // loop n_cnt
}
The console output of the program:
***************************************************************
* Running Example transient_ex2:
* mpirun -np 2 example-devel pipe-mesh.unv -pc_type bjacobi -sub_pc_type ilu -sub_pc_factor_levels 4 -sub_pc_factor_zeropivot 0 -ksp_right_pc -log_summary
***************************************************************
Running /workspace/libmesh/examples/transient/transient_ex2/.libs/lt-example-devel pipe-mesh.unv -pc_type bjacobi -sub_pc_type ilu -sub_pc_factor_levels 4 -sub_pc_factor_zeropivot 0 -ksp_right_pc -log_summary
Mesh file is: pipe-mesh.unv
*** Warning, This code is deprecated, and likely to be removed in future library versions! src/mesh/mesh_data.C, line 49, compiled Feb 5 2013 at 13:18:36 ***
Mesh Information:
mesh_dimension()=3
spatial_dimension()=3
n_nodes()=3977
n_local_nodes()=2037
n_elem()=3520
n_local_elem()=1760
n_active_elem()=3520
n_subdomains()=1
n_partitions()=2
n_processors()=2
n_threads()=1
processor_id()=0
EquationSystems
n_systems()=1
System #0, "Wave"
Type "Newmark"
Variables="p"
Finite Element Types="LAGRANGE", "JACOBI_20_00"
Infinite Element Mapping="CARTESIAN"
Approximation Orders="FIRST", "THIRD"
n_dofs()=3977
n_local_dofs()=2037
n_constrained_dofs()=0
n_local_constrained_dofs()=0
n_vectors()=9
n_matrices()=4
DofMap Sparsity
Average On-Processor Bandwidth <= 23.2034
Average Off-Processor Bandwidth <= 0.386724
Maximum On-Processor Bandwidth <= 27
Maximum Off-Processor Bandwidth <= 9
DofMap Constraints
Number of DoF Constraints = 0
Number of Node Constraints = 0
************************************************************************************************************************
*** WIDEN YOUR WINDOW TO 120 CHARACTERS. Use 'enscript -r -fCourier9' to print this document ***
************************************************************************************************************************
---------------------------------------------- PETSc Performance Summary: ----------------------------------------------
/workspace/libmesh/examples/transient/transient_ex2/.libs/lt-example-devel on a intel-12. named hbar.ices.utexas.edu with 2 processors, by benkirk Tue Feb 5 13:53:09 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.268e+01 1.00000 3.268e+01
Objects: 1.594e+03 1.00000 1.594e+03
Flops: 1.909e+10 1.06328 1.852e+10 3.704e+10
Flops/sec: 5.841e+08 1.06328 5.667e+08 1.133e+09
MPI Messages: 3.546e+03 1.00000 3.546e+03 7.093e+03
MPI Message Lengths: 2.881e+06 1.00000 8.123e+02 5.762e+06
MPI Reductions: 1.408e+04 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.2679e+01 100.0% 3.7039e+10 100.0% 7.093e+03 100.0% 8.123e+02 100.0% 1.408e+04 100.0%
------------------------------------------------------------------------------------------------------------------------
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
VecMDot 2023 1.0 4.9591e-02 2.2 3.24e+07 1.1 0.0e+00 0.0e+00 2.0e+03 0 0 0 0 14 0 0 0 0 14 1274
VecNorm 2623 1.0 2.3378e-02 1.9 1.07e+07 1.1 0.0e+00 0.0e+00 2.6e+03 0 0 0 0 19 0 0 0 0 19 892
VecScale 2623 1.0 3.4380e-03 1.2 5.34e+06 1.1 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 3034
VecCopy 1500 1.0 3.6411e-03 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
VecSet 3851 1.0 6.9311e-03 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 4200 1.0 2.3629e-02 1.1 1.59e+07 1.1 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 1313
VecMAXPY 2323 1.0 1.2114e-02 1.1 4.06e+07 1.1 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 6546
VecAssemblyBegin 1800 1.0 1.2624e-0110.6 0.00e+00 0.0 6.0e+02 5.8e+02 5.4e+03 0 0 8 6 38 0 0 8 6 38 0
VecAssemblyEnd 1800 1.0 1.1919e-03 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
VecScatterBegin 3223 1.0 9.8872e-03 1.1 0.00e+00 0.0 6.4e+03 8.1e+02 0.0e+00 0 0 91 91 0 0 0 91 91 0 0
VecScatterEnd 3223 1.0 2.0126e+00289.0 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00 3 0 0 0 0 3 0 0 0 0 0
VecNormalize 2323 1.0 2.4224e-02 1.6 1.42e+07 1.1 0.0e+00 0.0e+00 2.3e+03 0 0 0 0 16 0 0 0 0 17 1144
MatMult 2323 1.0 2.1346e+0016.2 2.17e+08 1.1 4.6e+03 7.8e+02 0.0e+00 3 1 66 63 0 3 1 66 63 0 198
MatMultAdd 600 1.0 3.7330e-02 1.1 5.72e+07 1.1 1.2e+03 7.8e+02 0.0e+00 0 0 17 16 0 0 0 17 16 0 2991
MatSolve 2623 1.0 9.5764e-01 1.0 2.28e+09 1.1 0.0e+00 0.0e+00 0.0e+00 3 12 0 0 0 3 12 0 0 0 4638
MatLUFactorNum 300 1.0 9.4695e+00 1.1 1.64e+10 1.1 0.0e+00 0.0e+00 0.0e+00 28 86 0 0 0 28 86 0 0 0 3364
MatILUFactorSym 300 1.0 2.0165e+01 1.1 0.00e+00 0.0 0.0e+00 0.0e+00 9.0e+02 60 0 0 0 6 60 0 0 0 6 0
MatAssemblyBegin 1207 1.0 1.7531e-02 1.4 0.00e+00 0.0 1.5e+01 1.1e+04 2.4e+03 0 0 0 3 17 0 0 0 3 17 0
MatAssemblyEnd 1207 1.0 5.4995e-02 1.0 0.00e+00 0.0 2.8e+01 2.0e+02 5.6e+01 0 0 0 0 0 0 0 0 0 0 0
MatGetRow 12222 1.1 1.3108e-03 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
MatGetRowIJ 300 1.0 9.2983e-05 1.7 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 300 1.0 7.1356e-03 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 6.0e+02 0 0 0 0 4 0 0 0 0 4 0
MatZeroEntries 10 1.0 5.3477e-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
MatAXPY 3 1.0 8.3830e-03 1.0 0.00e+00 0.0 1.2e+01 2.0e+02 4.8e+01 0 0 0 0 0 0 0 0 0 0 0
KSPGMRESOrthog 2023 1.0 6.1386e-02 1.8 6.47e+07 1.1 0.0e+00 0.0e+00 2.0e+03 0 0 0 0 14 0 0 0 0 14 2059
KSPSetUp 600 1.0 1.6849e-03 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 300 1.0 3.0923e+01 1.0 1.90e+10 1.1 4.6e+03 7.8e+02 6.1e+03 95100 66 63 44 95100 66 63 44 1193
PCSetUp 600 1.0 2.9737e+01 1.1 1.64e+10 1.1 0.0e+00 0.0e+00 1.5e+03 88 86 0 0 11 88 86 0 0 11 1071
PCSetUpOnBlocks 300 1.0 2.9736e+01 1.1 1.64e+10 1.1 0.0e+00 0.0e+00 1.5e+03 88 86 0 0 11 88 86 0 0 11 1071
PCApply 2623 1.0 9.8169e-01 1.0 2.28e+09 1.1 0.0e+00 0.0e+00 0.0e+00 3 12 0 0 0 3 12 0 0 0 4524
------------------------------------------------------------------------------------------------------------------------
Memory usage is given in bytes:
Object Type Creations Destructions Memory Descendants' Mem.
Reports information only for process 0.
--- Event Stage 0: Main Stage
Vector 343 343 5831192 0
Vector Scatter 8 8 8288 0
Index Set 916 916 3128944 0
IS L to G Mapping 1 1 564 0
Matrix 321 321 1575551704 0
Krylov Solver 2 2 19360 0
Preconditioner 2 2 1784 0
Viewer 1 0 0 0
========================================================================================================================
Average time to get PetscTime(): 9.53674e-08
Average time for MPI_Barrier(): 1.19209e-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:53:09 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=32.6895, Active time=32.4478 |
----------------------------------------------------------------------------------------------------------------
| 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() 1 0.0714 0.071352 0.0799 0.079859 0.22 0.25 |
| build_sparsity() 1 0.0675 0.067546 0.1557 0.155670 0.21 0.48 |
| create_dof_constraints() 1 0.0036 0.003630 0.0036 0.003630 0.01 0.01 |
| distribute_dofs() 1 0.0446 0.044560 0.1216 0.121571 0.14 0.37 |
| dof_indices() 10736 0.5158 0.000048 0.5158 0.000048 1.59 1.59 |
| prepare_send_list() 1 0.0002 0.000151 0.0002 0.000151 0.00 0.00 |
| reinit() 1 0.0761 0.076071 0.0761 0.076071 0.23 0.23 |
| |
| EquationSystems |
| build_solution_vector() 4 0.0286 0.007138 0.3654 0.091355 0.09 1.13 |
| |
| FE |
| compute_shape_functions() 1760 0.0406 0.000023 0.0406 0.000023 0.13 0.13 |
| init_shape_functions() 8 0.0004 0.000051 0.0004 0.000051 0.00 0.00 |
| |
| FEMap |
| compute_affine_map() 1440 0.0149 0.000010 0.0149 0.000010 0.05 0.05 |
| compute_map() 320 0.0074 0.000023 0.0074 0.000023 0.02 0.02 |
| init_reference_to_physical_map() 8 0.0004 0.000053 0.0004 0.000053 0.00 0.00 |
| |
| GMVIO |
| write_nodal_data() 4 0.1111 0.027764 0.1112 0.027812 0.34 0.34 |
| |
| Mesh |
| find_neighbors() 1 0.0642 0.064244 0.0644 0.064393 0.20 0.20 |
| read() 1 0.0006 0.000608 0.0434 0.043354 0.00 0.13 |
| renumber_nodes_and_elem() 2 0.0029 0.001426 0.0029 0.001426 0.01 0.01 |
| |
| MeshCommunication |
| broadcast() 1 0.0183 0.018305 0.0353 0.035289 0.06 0.11 |
| compute_hilbert_indices() 2 0.0330 0.016496 0.0330 0.016496 0.10 0.10 |
| find_global_indices() 2 0.0110 0.005511 0.0474 0.023692 0.03 0.15 |
| parallel_sort() 2 0.0027 0.001350 0.0028 0.001390 0.01 0.01 |
| |
| MeshOutput |
| write_equation_systems() 4 0.0003 0.000072 0.4772 0.119289 0.00 1.47 |
| |
| MetisPartitioner |
| partition() 1 0.0619 0.061857 0.0849 0.084877 0.19 0.26 |
| |
| NewmarkSystem |
| initial_conditions () 1 0.0000 0.000020 0.0000 0.000020 0.00 0.00 |
| update_rhs () 300 0.1049 0.000350 0.1049 0.000350 0.32 0.32 |
| update_u_v_a () 300 0.0172 0.000057 0.0172 0.000057 0.05 0.05 |
| |
| Parallel |
| allgather() 9 0.0004 0.000040 0.0004 0.000041 0.00 0.00 |
| broadcast() 6 0.0008 0.000129 0.0008 0.000129 0.00 0.00 |
| max(bool) 4 0.0000 0.000002 0.0000 0.000002 0.00 0.00 |
| max(scalar) 1689 0.0045 0.000003 0.0045 0.000003 0.01 0.01 |
| max(vector) 342 0.0022 0.000006 0.0046 0.000013 0.01 0.01 |
| min(bool) 2022 0.0043 0.000002 0.0043 0.000002 0.01 0.01 |
| min(scalar) 1682 0.0239 0.000014 0.0239 0.000014 0.07 0.07 |
| min(vector) 342 0.0027 0.000008 0.0054 0.000016 0.01 0.02 |
| probe() 12 0.0007 0.000062 0.0007 0.000062 0.00 0.00 |
| receive() 12 0.0001 0.000011 0.0009 0.000073 0.00 0.00 |
| send() 12 0.0001 0.000005 0.0001 0.000005 0.00 0.00 |
| send_receive() 16 0.0003 0.000021 0.0013 0.000082 0.00 0.00 |
| sum() 329 0.0089 0.000027 0.0116 0.000035 0.03 0.04 |
| |
| Parallel::Request |
| wait() 12 0.0000 0.000003 0.0000 0.000003 0.00 0.00 |
| |
| Partitioner |
| set_node_processor_ids() 1 0.0035 0.003549 0.0036 0.003627 0.01 0.01 |
| set_parent_processor_ids() 1 0.0027 0.002685 0.0027 0.002685 0.01 0.01 |
| |
| PetscLinearSolver |
| solve() 300 30.9668 0.103223 30.9668 0.103223 95.44 95.44 |
| |
| System |
| assemble() 1 0.0838 0.083778 0.2346 0.234590 0.26 0.72 |
| |
| UNVIO |
| count_elements() 1 0.0009 0.000854 0.0009 0.000854 0.00 0.00 |
| count_nodes() 1 0.0011 0.001139 0.0011 0.001139 0.00 0.00 |
| element_in() 1 0.0301 0.030145 0.0301 0.030145 0.09 0.09 |
| node_in() 1 0.0106 0.010607 0.0106 0.010607 0.03 0.03 |
----------------------------------------------------------------------------------------------------------------
| Totals: 21699 32.4478 100.00 |
----------------------------------------------------------------------------------------------------------------
***************************************************************
* Done Running Example transient_ex2:
* mpirun -np 2 example-devel pipe-mesh.unv -pc_type bjacobi -sub_pc_type ilu -sub_pc_factor_levels 4 -sub_pc_factor_zeropivot 0 -ksp_right_pc -log_summary
***************************************************************
Site Created By: libMesh Developers
Last modified: February 05 2013 19:53:09 UTC
Hosted By:
Transient Example 2 - The Newmark System and the Wave Equation
This is the eighth example program. It builds on the previous example programs. It introduces the NewmarkSystem class. In this example the wave equation is solved using the time integration scheme provided by the NewmarkSystem class.
This example comes with a cylindrical mesh given in the universal file pipe-mesh.unv. The mesh contains HEX8 and PRISM6 elements.
C++ include files that we need