// @sect3{The <code>FluidStructureProblem</code> implementation}
+ // @sect4{Constructors and helper functions}
+
// Let's now get to the implementation of the
// primary class of this program. The first
- // few functions are the constructor and
-
-
+ // few functions are the constructor and the
+ // helper functions that can be used to
+ // determine which part of the domain a cell
+ // is in. Given the discussion of these
+ // topics in the introduction, their
+ // implementation is rather obvious. In the
+ // constructor, note that we have to
+ // construct the hp::FECollection object from
+ // the base elements for Stokes and
+ // elasticity; using the
+ // hp::FECollection::push_back function
+ // assigns them spots zero and one in this
+ // collection, an order that we have to
+ // remember and use consistently in the rest
+ // of the program.
template <int dim>
FluidStructureProblem<dim>::
FluidStructureProblem (const unsigned int stokes_degree,
}
-
+ // @sect4{Meshes and assigning subdomains}
+
+ // The next pair of functions deals with
+ // generating a mesh and making sure all
+ // flags that denote subdomains are
+ // correct. <code>make_grid</code>, as
+ // discussed in the introduction, generates
+ // an $8\times 8$ mesh (or an $8\times
+ // 8\times 8$ mesh in 3d) to make sure that
+ // each coarse mesh cell is completely within
+ // one of the subdomains. After generating
+ // this mesh, we loop over its boundary and
+ // set the boundary indicator to one at the
+ // top boundary, the only place where we set
+ // nonzero Dirichlet boundary
+ // conditions. After this, we loop again over
+ // all cells to set the material indicator
+ // — used to denote which part of the
+ // domain we are in, to either the fluid or
+ // solid indicator.
template <int dim>
void
FluidStructureProblem<dim>::make_grid ()
{
-// not quite what we want...
GridGenerator::subdivided_hyper_cube (triangulation, 8, -1, 1);
+
for (typename Triangulation<dim>::active_cell_iterator
cell = triangulation.begin_active();
cell != triangulation.end(); ++cell)
}
-
+ // The second part of this pair of functions
+ // determines which finite element to use on
+ // each cell. Above we have set the material
+ // indicator for each coarse mesh cell, and
+ // as mentioned in the introduction, this
+ // information is inherited from mother to
+ // child cell upon mesh refinement.
+ //
+ // In other words, whenever we have refined
+ // (or created) the mesh, we can rely on the
+ // material indicators to be a correct
+ // description of which part of the domain a
+ // cell is in. We then use this to set the
+ // active FE index of the cell to the
+ // corresponding element of the
+ // hp::FECollection member variable of this
+ // class: zero for fluid cells, one for solid
+ // cells.
template <int dim>
void
FluidStructureProblem<dim>::set_active_fe_indices ()
for (typename hp::DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active();
cell != dof_handler.end(); ++cell)
- if (cell_is_in_fluid_domain(cell))
- cell->set_active_fe_index (0);
- else if (cell_is_in_solid_domain(cell))
- cell->set_active_fe_index (1);
- else
- Assert (false, ExcNotImplemented());
+ {
+ if (cell_is_in_fluid_domain(cell))
+ cell->set_active_fe_index (0);
+ else if (cell_is_in_solid_domain(cell))
+ cell->set_active_fe_index (1);
+ else
+ Assert (false, ExcNotImplemented());
+ }
}
-
+ // @sect4{<code>FluidStructureProblem::setup_dofs</code>}
+
+ // The next step is to setup the data
+ // structures for the linear system. To this
+ // end, we first have to set the active FE
+ // indices with the function immediately
+ // above, then distribute degrees of freedom,
+ // and then determine constraints on the
+ // linear system. The latter includes hanging
+ // node constraints as usual, but also the
+ // inhomogenous boundary values at the top
+ // fluid boundary, and zero boundary values
+ // along the perimeter of the solid
+ // subdomain.
template <int dim>
-void FluidStructureProblem<dim>::setup_dofs ()
+void
+FluidStructureProblem<dim>::setup_dofs ()
{
set_active_fe_indices ();
dof_handler.distribute_dofs (fe_collection);
StokesBoundaryValues<dim>(),
constraints,
velocity_mask);
+
std::vector<bool> elasticity_mask (dim+1+dim, false);
for (unsigned int d=dim+1; d<dim+1+dim; ++d)
elasticity_mask[d] = true;
elasticity_mask);
}
- // make sure velocity is zero at
- // the interface
+ // There are more constraints we have to
+ // handle, though: we have to make sure
+ // that the velocity is zero at the
+ // interface between fluid and solid. The
+ // following piece of code was already
+ // presented in the introduction:
{
std::vector<unsigned int> local_face_dof_indices (stokes_fe.dofs_per_face);
for (typename hp::DoFHandler<dim>::active_cell_iterator
face_is_on_interface = true;
else if (cell->neighbor(f)->has_children() == true)
{
- // neighbor does
- // have
- // children. see if
- // any of the cells
- // on the other
- // side are elastic
for (unsigned int sf=0; sf<cell->face(f)->n_children(); ++sf)
- if (cell_is_in_solid_domain (cell->neighbor_child_on_subface(f, sf)))
+ if (cell_is_in_solid_domain (cell->neighbor_child_on_subface
+ (f, sf)))
{
face_is_on_interface = true;
break;
}
}
-
+ // At the end of all this, we can declare
+ // to the constraints object that we now
+ // have all constraints ready to go and
+ // that the object can rebuild its internal
+ // data structures for better efficiency:
constraints.close ();
std::cout << " Number of active cells: "
<< dof_handler.n_dofs()
<< std::endl;
+ // The rest of this function is standard:
+ // Create a sparsity pattern and use it to
+ // initialize the matrix; then also set
+ // vectors to their correct sizes.
{
CompressedSimpleSparsityPattern csp (dof_handler.n_dofs(),
dof_handler.n_dofs());
+ // @sect4{<code>FluidStructureProblem::assemble_system</code>}
+
+ // Following is the central function of this
+ // program: the one that assembles the linear
+ // system. It has a long section of setting
+ // up auxiliary functions at the beginning:
+ // from creating the quadrature formulas and
+ // setting up the FEValues, FEFaceValues and
+ // FESubfaceValues objects necessary to
+ // integrate the cell terms as well as the
+ // interface terms for the case where cells
+ // along the interface come together at same
+ // size or with differing levels of
+ // refinement...
template <int dim>
void FluidStructureProblem<dim>::assemble_system ()
{
update_JxW_values |
update_gradients);
- const QGauss<dim-1> face_quadrature(std::max (stokes_degree+2,
- elasticity_degree+2));
+ const QGauss<dim-1> common_face_quadrature(std::max (stokes_degree+2,
+ elasticity_degree+2));
FEFaceValues<dim> stokes_fe_face_values (stokes_fe,
- face_quadrature,
+ common_face_quadrature,
update_JxW_values |
update_normal_vectors |
update_gradients);
FEFaceValues<dim> elasticity_fe_face_values (elasticity_fe,
- face_quadrature,
+ common_face_quadrature,
update_values);
FESubfaceValues<dim> stokes_fe_subface_values (stokes_fe,
- face_quadrature,
+ common_face_quadrature,
update_JxW_values |
update_normal_vectors |
update_gradients);
FESubfaceValues<dim> elasticity_fe_subface_values (elasticity_fe,
- face_quadrature,
+ common_face_quadrature,
update_values);
- const unsigned int stokes_dofs_per_cell = stokes_fe.dofs_per_cell;
- const unsigned int elasticity_dofs_per_cell = elasticity_fe.dofs_per_cell;
+ // ...to objects that are needed to
+ // describe the local contributions to the
+ // global linear system...
+ const unsigned int stokes_dofs_per_cell = stokes_fe.dofs_per_cell;
+ const unsigned int elasticity_dofs_per_cell = elasticity_fe.dofs_per_cell;
- FullMatrix<double> local_matrix;
- FullMatrix<double> local_interface_matrix (elasticity_dofs_per_cell,
- stokes_dofs_per_cell);
- Vector<double> local_rhs;
+ FullMatrix<double> local_matrix;
+ FullMatrix<double> local_interface_matrix (elasticity_dofs_per_cell,
+ stokes_dofs_per_cell);
+ Vector<double> local_rhs;
std::vector<unsigned int> local_dof_indices;
std::vector<unsigned int> neighbor_dof_indices (stokes_dofs_per_cell);
- const RightHandSide<dim> right_hand_side;
+ const RightHandSide<dim> right_hand_side;
- const FEValuesExtractors::Vector velocities (0);
- const FEValuesExtractors::Scalar pressure (dim);
- const FEValuesExtractors::Vector displacements (dim+1);
+ // ...to variables that allow us to extract
+ // certain components of the shape
+ // functions and cache their values rather
+ // than having to recompute them at every
+ // quadrature point:
+ const FEValuesExtractors::Vector velocities (0);
+ const FEValuesExtractors::Scalar pressure (dim);
+ const FEValuesExtractors::Vector displacements (dim+1);
std::vector<SymmetricTensor<2,dim> > stokes_phi_grads_u (stokes_dofs_per_cell);
std::vector<double> stokes_div_phi_u (stokes_dofs_per_cell);
std::vector<Tensor<2,dim> > elasticity_phi_grad (elasticity_dofs_per_cell);
std::vector<double> elasticity_phi_div (elasticity_dofs_per_cell);
- std::vector<Tensor<1,dim> > elasticity_phi (elasticity_dofs_per_cell);
-
+ std::vector<Tensor<1,dim> > elasticity_phi (elasticity_dofs_per_cell);
+
+ // Then comes the main loop over all cells
+ // and, as in step-27, the initialization
+ // of the hp::FEValues object for the
+ // current cell and the extraction of a
+ // FEValues object that is appropriate for
+ // the current cell:
typename hp::DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
cell->get_fe().dofs_per_cell);
local_rhs.reinit (cell->get_fe().dofs_per_cell);
+ // With all of this done, we continue
+ // to assemble the cell terms for cells
+ // that are part of the Stokes and
+ // elastic regions. While we could in
+ // principle do this in one formula, in
+ // effect implementing the one bilinear
+ // form stated in the introduction, we
+ // realize that our finite element
+ // spaces are chosen in such a way that
+ // on each cell, one set of variables
+ // (either velocities and pressure, or
+ // displacements) are always zero, and
+ // consequently a more efficient way of
+ // computing local integrals is to do
+ // only what's necessary based on an
+ // <code>if</code> clause that tests
+ // which part of the domain we are in.
+ //
+ // The actual computation of the local
+ // matrix is the same as in step-22 as
+ // well as that given in the @ref
+ // vector_valued documentation module
+ // for the elasticity equations:
if (cell_is_in_fluid_domain (cell))
{
const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
}
}
+ // Once we have the contributions from
+ // cell integrals, we copy them into
+ // the global matrix (taking care of
+ // constraints right away, through the
+ // ConstraintMatrix::distribute_local_to_global
+ // function). Note that we have not
+ // written anything into the
+ // <code>local_rhs</code> variable,
+ // though we still need to pass it
+ // along since the elimination of
+ // nonzero boundary values requires the
+ // modification of local and
+ // consequently also global right hand
+ // side values:
local_dof_indices.resize (cell->get_fe().dofs_per_cell);
cell->get_dof_indices (local_dof_indices);
-
- // local_rhs==0, but need to do
- // this here because of
- // boundary values
constraints.distribute_local_to_global (local_matrix, local_rhs,
local_dof_indices,
system_matrix, system_rhs);
- // see about face terms
+ // The more interesting part of this
+ // function is where we see about face
+ // terms along the interface between
+ // the two subdomains. To this end, we
+ // first have to make sure that we only
+ // assemble them once even though a
+ // loop over all faces of all cells
+ // would encounter each part of the
+ // interface twice. We arbitrarily make
+ // the decision that we will only
+ // evaluate interface terms if the
+ // current cell is part of the solid
+ // subdomain and if, consequently, a
+ // face is not at the boundary and the
+ // potential neighbor behind it is part
+ // of the fluid domain. Let's start
+ // with these conditions:
if (cell_is_in_solid_domain (cell))
- // we are on a solid cell
for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
if (cell->at_boundary(f) == false)
{
+ // At this point we know that
+ // the current cell is a
+ // candidate for integration
+ // and that a neighbor behind
+ // face <code>f</code>
+ // exists. There are now three
+ // possibilities:
+ //
+ // - The neighbor is at the
+ // same refinement level and
+ // has no children.
+ // - The neighbor has children.
+ // - The neighbor is coarser.
+ //
+ // In all three cases, we are
+ // only interested in it if it
+ // is part of the fluid
+ // subdomain. So let us start
+ // with the first and simplest
+ // case: if the neighbor is at
+ // the same level, has no
+ // children, and is a fluid
+ // cell, then the two cells
+ // share a boundary that is
+ // part of the interface along
+ // which we want to integrate
+ // interface terms. All we have
+ // to do is initialize two
+ // FEFaceValues object with the
+ // current face and the face of
+ // the neighboring cell (note
+ // how we find out which face
+ // of the neighboring cell
+ // borders on the current cell)
+ // and pass things off to the
+ // function that evaluates the
+ // interface terms (the third
+ // through fifth arguments to
+ // this function provide it
+ // with scratch arrays). The
+ // result is then again copied
+ // into the global matrix,
+ // using a function that knows
+ // that the DoF indices of rows
+ // and columns of the local
+ // matrix result from different
+ // cells:
if ((cell->neighbor(f)->level() == cell->level())
&&
(cell->neighbor(f)->has_children() == false)
&&
cell_is_in_fluid_domain (cell->neighbor(f)))
{
- // same size
- // neighbors;
- // neighbor is
- // fluid cell
elasticity_fe_face_values.reinit (cell, f);
stokes_fe_face_values.reinit (cell->neighbor(f),
cell->neighbor_of_neighbor(f));
neighbor_dof_indices,
system_matrix);
}
+
+ // The second case is if the
+ // neighbor has further
+ // children. In that case, we
+ // have to loop over all the
+ // children of the neighbor to
+ // see if they are part of the
+ // fluid subdomain. If they
+ // are, then we integrate over
+ // the common interface, which
+ // is a face for the neighbor
+ // and a subface of the current
+ // cell, requiring us to use an
+ // FEFaceValues for the
+ // neighbor and an
+ // FESubfaceValues for the
+ // current cell:
else if ((cell->neighbor(f)->level() == cell->level())
&&
(cell->neighbor(f)->has_children() == true))
{
- // neighbor has children. loop over
- // the cells adjacent to the commone
- // interface and see which subdomain
- // they belong to
- for (unsigned int subface=0; subface<cell->face(f)->n_children(); ++subface)
- if (cell_is_in_fluid_domain (cell->neighbor_child_on_subface (f, subface)))
+ for (unsigned int subface=0;
+ subface<cell->face(f)->n_children();
+ ++subface)
+ if (cell_is_in_fluid_domain (cell->neighbor_child_on_subface
+ (f, subface)))
{
elasticity_fe_subface_values.reinit (cell,
f,
stokes_fe_face_values.reinit (cell->neighbor_child_on_subface (f, subface),
cell->neighbor_of_neighbor(f));
- assemble_interface_term (elasticity_fe_subface_values, stokes_fe_face_values,
- elasticity_phi, stokes_phi_grads_u, stokes_phi_p,
+ assemble_interface_term (elasticity_fe_subface_values,
+ stokes_fe_face_values,
+ elasticity_phi,
+ stokes_phi_grads_u, stokes_phi_p,
local_interface_matrix);
- cell->neighbor_child_on_subface (f, subface)->get_dof_indices (neighbor_dof_indices);
+ cell->neighbor_child_on_subface (f, subface)
+ ->get_dof_indices (neighbor_dof_indices);
constraints.distribute_local_to_global(local_interface_matrix,
local_dof_indices,
neighbor_dof_indices,
system_matrix);
}
}
+
+ // The last option is that the
+ // neighbor is coarser. In that
+ // case we have to use an
+ // FESubfaceValues object for
+ // the neighbor and a
+ // FEFaceValues for the current
+ // cell; the rest is the same
+ // as before:
else if (cell->neighbor_is_coarser(f)
&&
cell_is_in_fluid_domain(cell->neighbor(f)))
{
- // neighbor is coarser
elasticity_fe_face_values.reinit (cell, f);
stokes_fe_subface_values.reinit (cell->neighbor(f),
cell->neighbor_of_coarser_neighbor(f).first,
cell->neighbor_of_coarser_neighbor(f).second);
- assemble_interface_term (elasticity_fe_face_values, stokes_fe_subface_values,
- elasticity_phi, stokes_phi_grads_u, stokes_phi_p,
+ assemble_interface_term (elasticity_fe_face_values,
+ stokes_fe_subface_values,
+ elasticity_phi,
+ stokes_phi_grads_u, stokes_phi_p,
local_interface_matrix);
cell->neighbor(f)->get_dof_indices (neighbor_dof_indices);
+ // In the function that assembles the global
+ // system, we passed computing interface
+ // terms to a separate function we discuss
+ // here. The key is that even though we can't
+ // predict the combination of FEFaceValues
+ // and FESubfaceValues objects, they are both
+ // derived from the FEFaceValuesBase class
+ // and consequently we don't have to care:
+ // the function is simply called with two
+ // such objects denoting the values of the
+ // shape functions on the quadrature points
+ // of the two sides of the face. We then do
+ // what we always do: we fill the scratch
+ // arrays with the values of shape functions
+ // and their derivatives, and then loop over
+ // all entries of the matrix to compute the
+ // local integrals. The details of the
+ // bilinear form we evaluate here are given
+ // in the introduction.
template <int dim>
void
-FluidStructureProblem<dim>::assemble_interface_term (const FEFaceValuesBase<dim> &elasticity_fe_face_values,
- const FEFaceValuesBase<dim> &stokes_fe_face_values,
- std::vector<Tensor<1,dim> > &elasticity_phi,
- std::vector<SymmetricTensor<2,dim> > &stokes_phi_grads_u,
- std::vector<double> &stokes_phi_p,
- FullMatrix<double> &local_interface_matrix) const
+FluidStructureProblem<dim>::
+assemble_interface_term (const FEFaceValuesBase<dim> &elasticity_fe_face_values,
+ const FEFaceValuesBase<dim> &stokes_fe_face_values,
+ std::vector<Tensor<1,dim> > &elasticity_phi,
+ std::vector<SymmetricTensor<2,dim> > &stokes_phi_grads_u,
+ std::vector<double> &stokes_phi_p,
+ FullMatrix<double> &local_interface_matrix) const
{
Assert (stokes_fe_face_values.n_quadrature_points ==
elasticity_fe_face_values.n_quadrature_points,
ExcInternalError());
-
+ const unsigned int n_face_quadrature_points
+ = elasticity_fe_face_values.n_quadrature_points;
+
const FEValuesExtractors::Vector velocities (0);
const FEValuesExtractors::Scalar pressure (dim);
const FEValuesExtractors::Vector displacements (dim+1);
local_interface_matrix = 0;
- for (unsigned int q=0; q<elasticity_fe_face_values.n_quadrature_points; ++q)
+ for (unsigned int q=0; q<n_face_quadrature_points; ++q)
{
const Tensor<1,dim> normal_vector = stokes_fe_face_values.normal_vector(q);
}
+ // @sect4{<code>FluidStructureProblem::solve</code>}
+
+ // As discussed in the introduction, we use a
+ // rather trivial solver here: we just pass
+ // the linear system off to the
+ // SparseDirectUMFPACK direct solver (see,
+ // for example, step-29). The only thing we
+ // have to do after solving is ensure that
+ // hanging node and boundary value
+ // constraints are correct.
template <int dim>
void
FluidStructureProblem<dim>::solve ()
+ // @sect4{<code>FluidStructureProblem::output_results</code>}
+ // Generating graphical output is rather
+ // trivial here: all we have to do is
+ // identify which components of the solution
+ // vector belong to scalars and/or vectors
+ // (see, for example, step-22 for a previous
+ // example), and then pass it all on to the
+ // DataOut class (with the second template
+ // argument equal to hp::DoFHandler instead
+ // of the usual default DoFHandler):
template <int dim>
void
-FluidStructureProblem<dim>::output_results (const unsigned int refinement_cycle) const
+FluidStructureProblem<dim>::
+output_results (const unsigned int refinement_cycle) const
{
std::vector<std::string> solution_names (dim, "velocity");
solution_names.push_back ("pressure");
}
-
+ // @sect4{<code>FluidStructureProblem::refine_mesh</code>}
+
+ // The next step is to refine the mesh. As
+ // was discussed in the introduction, this is
+ // a bit tricky primarily because the fluid
+ // and the solid subdomains use variables
+ // that have different physical dimensions
+ // and for which the absolute magnitude of
+ // error estimates is consequently not
+ // directly comparable. We will therefore
+ // have to scale them. At the top of the
+ // function, we therefore first compute error
+ // estimates for the different variables
+ // separately (using the velocities but not
+ // the pressure for the fluid domain, and the
+ // displacements in the solid domain):
template <int dim>
void
FluidStructureProblem<dim>::refine_mesh ()
{
- Vector<float> stokes_estimated_error_per_cell (triangulation.n_active_cells());
- Vector<float> elasticity_estimated_error_per_cell (triangulation.n_active_cells());
- Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
+ Vector<float>
+ stokes_estimated_error_per_cell (triangulation.n_active_cells());
+ Vector<float>
+ elasticity_estimated_error_per_cell (triangulation.n_active_cells());
- const QGauss<dim-1> stokes_quadrature(stokes_degree+2);
- const QGauss<dim-1> elasticity_quadrature(elasticity_degree+2);
+ const QGauss<dim-1> stokes_face_quadrature(stokes_degree+2);
+ const QGauss<dim-1> elasticity_face_quadrature(elasticity_degree+2);
hp::QCollection<dim-1> face_q_collection;
- face_q_collection.push_back (stokes_quadrature);
- face_q_collection.push_back (elasticity_quadrature);
+ face_q_collection.push_back (stokes_face_quadrature);
+ face_q_collection.push_back (elasticity_face_quadrature);
std::vector<bool> stokes_component_mask (dim+1+dim, false);
for (unsigned int d=0; d<dim; ++d)
elasticity_estimated_error_per_cell,
elasticity_component_mask);
- stokes_estimated_error_per_cell /= 0.25 * stokes_estimated_error_per_cell.l2_norm();
- elasticity_estimated_error_per_cell /= elasticity_estimated_error_per_cell.l2_norm();
+ // We then normalize error estimates by
+ // dividing by their norm and scale the
+ // fluid error indicators by a factor of 4
+ // as discussed in the introduction. The
+ // results are then added together into a
+ // vector that contains error indicators
+ // for all cells:
+ stokes_estimated_error_per_cell
+ *= 4 . / stokes_estimated_error_per_cell.l2_norm();
+ elasticity_estimated_error_per_cell
+ *= 1. / elasticity_estimated_error_per_cell.l2_norm();
+
+ Vector<float>
+ estimated_error_per_cell (triangulation.n_active_cells());
+
estimated_error_per_cell += stokes_estimated_error_per_cell;
estimated_error_per_cell += elasticity_estimated_error_per_cell;
+ // The second to last part of the function,
+ // before actually refining the mesh,
+ // involves a heuristic that we have
+ // already mentioned in the introduction:
+ // because the solution is discontinuous,
+ // the KellyErrorEstimator class gets all
+ // confused about cells that sit at the
+ // boundary between subdomains: it believes
+ // that the error is large there because
+ // the jump in the gradient is large, even
+ // though this is entirely expected and a
+ // feature that is in fact present in the
+ // exact solution as well and therefore not
+ // indicative of any numerical error.
+ //
+ // Consequently, we set the error
+ // indicators to zero for all cells at the
+ // interface; the conditions determining
+ // which cells this affects are slightly
+ // awkward because we have to account for
+ // the possibility of adaptively refined
+ // meshes, meaning that the neighboring
+ // cell can be coarser than the current
+ // one, or could in fact be refined some
+ // more. The structure of these nested
+ // conditions is much the same as we
+ // encountered when assembling interface
+ // terms in <code>assemble_system</code>.
{
unsigned int cell_index = 0;
for (typename hp::DoFHandler<dim>::active_cell_iterator
+ // @sect4{<code>FluidStructureProblem::run</code>}
+
+ // This is, as usual, the function that
+ // controls the overall flow of operation. If
+ // you've read through tutorial programs
+ // step-1 through step-6, for example, then
+ // you are already quite familiar with the
+ // following structure:
template <int dim>
void FluidStructureProblem<dim>::run ()
{
+ // @sect4{The <code>main()</code> function}
+
+ // This, final, function contains pretty much
+ // exactly what most of the other tutorial
+ // programs have:
int main ()
{
try
projects / dealii.git / commitdiff