#include <iostream>
#include <fstream>
#include <algorithm>
#include <stdio.h>
#include <math.h>
Basic include file needed for the mesh functionality.
#include "libmesh.h"
#include "mesh.h"
#include "gmv_io.h"
#include "newmark_system.h"
#include "equation_systems.h"
Define the Finite Element object.
#include "fe.h"
Define Gauss quadrature rules.
#include "quadrature_gauss.h"
Define useful datatypes for finite element
#include "dense_matrix.h"
#include "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 "sparse_matrix.h"
#include "numeric_vector.h"
Define the DofMap, which handles degree of freedom
indexing.
#include "dof_map.h"
The definition of a vertex associated with a Mesh.
#include "node.h"
The definition of a geometric element
#include "elem.h"
Defines the MeshData class, which allows you to store
data about the mesh when reading in files, etc.
#include "mesh_data.h"
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 Petsc, like in example 2.
libMesh::init (argc, argv);
Braces are used to force object scope.
{
Check for proper usage.
if (argc < 2)
{
std::cerr << "Usage: " << argv[0] << " [meshfile]"
<< std::endl;
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). Print a warning to the user if PETSc
is not available, or if they are using LasPack solvers.
#ifdef HAVE_PETSC
if ((libMesh::on_command_line("--use-laspack")) ||
(libMesh::on_command_line("--disable-petsc")))
#endif
{
std::cerr << "WARNING! It appears you are using the\n"
<< "LasPack solvers. ex8 is known not to converge\n"
<< "using LasPack, but should work OK with PETSc.\n"
<< "If possible, download and install the PETSc\n"
<< "library from www-unix.mcs.anl.gov/petsc/petsc-2/\n"
<< std::endl;
}
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;
For now, restrict to dim=3, though this
may easily be changed, see example 4
const unsigned int dim = 3;
Create a dim-dimensional mesh.
Mesh mesh (dim);
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.
{
Creates 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");
Adds 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.;
Store the current time as an
\p EquationSystems parameter, so that
\p fill_dirichlet_bc() can access it.
equation_systems.parameters.set<Real>("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);
Use a handy reference to this system
NewmarkSystem& t_system = equation_systems.get_system<NewmarkSystem> ("Wave");
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
Real t_time = 0.;
res_out << "# pressure at node " << res_node_no << "\n"
<< "# time\tpressure\n"
<< t_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_time += delta_t;
equation_systems.parameters.set<Real>("time") = t_time;
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];
sprintf (buf, "out.%03d.gmv", time_step);
GMVIO(mesh).write_equation_systems (buf,
equation_systems);
}
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_time << "\t"
<< global_displacement[dof_no]
<< std::endl;
}
}
All done.
return libMesh::close ();
}
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.
assert (system_name == "Wave");
Get a constant reference to the mesh object.
const Mesh& 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<unsigned int> dof_indices;
Now we will loop over all the elements in the mesh.
We will compute the element matrix and right-hand-side
contribution.
const_elem_iterator el (mesh.elements_begin());
const const_elem_iterator end_el (mesh.elements_end());
MeshBase::const_element_iterator el = mesh.elements_begin();
const MeshBase::const_element_iterator end_el = mesh.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
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.
assert (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 Mesh& 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");
Ge the current time from \p EquationSystems
const Real current_time = es.parameters.get<Real>("time");
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 (current_time < .002 )
p_value = sin(2*pi*current_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 program without comments:
#include <iostream>
#include <fstream>
#include <algorithm>
#include <stdio.h>
#include <math.h>
#include "libmesh.h"
#include "mesh.h"
#include "gmv_io.h"
#include "newmark_system.h"
#include "equation_systems.h"
#include "fe.h"
#include "quadrature_gauss.h"
#include "dense_matrix.h"
#include "dense_vector.h"
#include "sparse_matrix.h"
#include "numeric_vector.h"
#include "dof_map.h"
#include "node.h"
#include "elem.h"
#include "mesh_data.h"
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)
{
libMesh::init (argc, argv);
{
if (argc < 2)
{
std::cerr << "Usage: " << argv[0] << " [meshfile]"
<< std::endl;
error();
}
else
{
std::cout << "Running " << argv[0];
for (int i=1; i<argc; i++)
std::cout << " " << argv[i];
std::cout << std::endl << std::endl;
}
#ifdef HAVE_PETSC
if ((libMesh::on_command_line("--use-laspack")) ||
(libMesh::on_command_line("--disable-petsc")))
#endif
{
std::cerr << "WARNING! It appears you are using the\n"
<< "LasPack solvers. ex8 is known not to converge\n"
<< "using LasPack, but should work OK with PETSc.\n"
<< "If possible, download and install the PETSc\n"
<< "library from www-unix.mcs.anl.gov/petsc/petsc-2/\n"
<< std::endl;
}
std::string mesh_file = argv[1];
std::cout << "Mesh file is: " << mesh_file << std::endl;
const unsigned int dim = 3;
Mesh mesh (dim);
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.;
equation_systems.parameters.set<Real>("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);
NewmarkSystem& t_system = equation_systems.get_system<NewmarkSystem> ("Wave");
t_system.assemble();
Real t_time = 0.;
res_out << "# pressure at node " << res_node_no << "\n"
<< "# time\tpressure\n"
<< t_time << "\t" << 0 << std::endl;
for (unsigned int time_step=0; time_step<n_time_steps; time_step++)
{
t_time += delta_t;
equation_systems.parameters.set<Real>("time") = t_time;
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];
sprintf (buf, "out.%03d.gmv", time_step);
GMVIO(mesh).write_equation_systems (buf,
equation_systems);
}
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_time << "\t"
<< global_displacement[dof_no]
<< std::endl;
}
}
return libMesh::close ();
}
void assemble_wave(EquationSystems& es,
const std::string& system_name)
{
assert (system_name == "Wave");
const Mesh& 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<unsigned int> dof_indices;
MeshBase::const_element_iterator el = mesh.elements_begin();
const MeshBase::const_element_iterator end_el = mesh.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)
{
assert (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 Mesh& mesh = es.get_mesh();
const Real pi = libMesh::pi;
const bool do_for_matrix =
es.parameters.get<bool>("Newmark set BC for Matrix");
const Real current_time = es.parameters.get<Real>("time");
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 (current_time < .002 )
p_value = sin(2*pi*current_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 ./ex8-devel pipe-mesh.unv
***************************************************************
Running ./ex8-devel pipe-mesh.unv
Mesh file is: pipe-mesh.unv
Mesh Information:
mesh_dimension()=3
spatial_dimension()=3
n_nodes()=3977
n_elem()=3520
n_local_elem()=3520
n_active_elem()=3520
n_subdomains()=1
n_processors()=1
processor_id()=0
EquationSystems
n_systems()=1
System "Wave"
Type "Newmark"
Variables="p"
Finite Element Types="LAGRANGE"
Approximation Orders="FIRST"
n_dofs()=3977
n_local_dofs()=3977
n_constrained_dofs()=0
n_vectors()=9
***************************************************************
* Done Running Example ./ex8-devel pipe-mesh.unv
***************************************************************
Site Created By: libMesh Developers
Last modified: October 22 2008 00:23:47.
Hosted By:
Example 8 - 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