Table of contents | |
---|---|
This program was contributed by Markus Bürg.
In this example we consider how to use periodic boundary conditions in deal.II. Periodic boundary conditions are a typical approach when one wants to solve some equation on a representative piece of a larger domain that repeats in one or more direction; an example is the simulation of the electronic structure of photonic crystals, because they have a lattice-like structure and, thus, it often suffices to do the actual computation on only one cell of the lattice. To be able to proceed this way one has to assume that the computation can be periodically extended to the other cells. This requires the solution to be periodic with respect to the cells. Hence the solution has to obtain the same nodal values on opposite parts of the boundary. In the figure below we show this concept in two space-dimensions. There, all dashed faces with the same color should have the same boundary values:
To keep things simple, in this tutorial we will consider an academic, simplified problem that allows us to focus on only that part that we are interested in here, namely how to set up periodic boundary conditions. Specifically, we solve the Poisson problem on a domain, where the left and right parts of the boundary are identified. Let and consider the problem
Note that the source term is not symmetric and so the solution would not be periodic unless this is explicitly enforced.
The way one has to see these periodic boundary conditions is as follows: Assume for a moment (as we do in this program) that we have a uniformly refined mesh. Then, after discretization there are a number of nodes (degrees of freedom) with indices
on the left boundary of the domain, and a second set of nodes at the right boundary
. Since we have assumed that the mesh is uniformly refined, there is exactly one node
for each
so that
, i.e. the two of them match with respect to the periodicity. We will then write that
(and, if you want,
). If now
are the unknowns of our discretized problem, then the periodic boundary condition boils down to the following set of constraints:
Now, this is exactly the sort of constraint that the ConstraintMatrix class, first introduced in step-6, handles and can enforce in a linear system. Consequently, the main point of this program is how we fill the ConstraintMatrix object that stores these constraints, and how this is applied to the resulting linear system.
The code for solving this problem is simple and based on step-3 since we want to focus on the implementation of the periodic boundary conditions. The code could be much more sophisticated, of course. For example, we could want to enforce periodic boundary conditions for adaptively refined meshes in which there is no longer a one-to-one relationship between degrees of freedom. We will discuss this at the end of the results section of this program.
The include files are already known. The one critical for the current program is the one that contains the ConstraintMatrix in the lac/
directory:
#include <base/function.h> #include <base/quadrature_lib.h> #include <lac/constraint_matrix.h> #include <lac/precondition.h> #include <lac/solver_cg.h> #include <lac/solver_control.h> #include <lac/sparse_matrix.h> #include <lac/sparsity_pattern.h> #include <lac/compressed_sparsity_pattern.h> #include <grid/grid_generator.h> #include <grid/tria.h> #include <dofs/dof_accessor.h> #include <dofs/dof_handler.h> #include <dofs/dof_tools.h> #include <fe/fe_q.h> #include <fe/fe_values.h> #include <numerics/data_out.h> #include <numerics/vectors.h> #include <fstream> using namespace dealii;
LaplaceProblem
classThe class LaplaceProblem
is the main class of this problem. As mentioned in the introduction, it is fashioned after the corresponding class in step-3. Correspondingly, the documentation from that tutorial program applies here as well. The only new member variable is the constraints
variables that will hold the constraints from the periodic boundary condition. We will initialize it in the make_periodicity_constraints()
function which we call from make_grid_and_dofs()
.
class LaplaceProblem { public: LaplaceProblem (); void run (); private: Triangulation<2> triangulation; FE_Q<2> fe; DoFHandler<2> dof_handler; ConstraintMatrix constraints; SparsityPattern sparsity_pattern; SparseMatrix<double> system_matrix; Vector<double> system_rhs; Vector<double> solution; void assemble_system (); void output_results (); void make_grid_and_dofs (); void make_periodicity_constraints (); void solve (); };
RightHandSide
classThe following implements the right hand side function discussed in the introduction. Its implementation is obvious given what has been shown in step-4 before:
class RightHandSide: public Function<2> { public: RightHandSide (); virtual double value (const Point<2>& p, const unsigned int component = 0) const; }; RightHandSide::RightHandSide () : Function<2> () {} double RightHandSide::value (const Point<2>&p, const unsigned int) const { return (std::cos (2 * numbers::PI * p(0)) * std::exp (- 2 * p(0)) * std::cos (2 * numbers::PI * p(1)) * std::exp (- 2 * p(1))); }
LaplaceProblem
classThe first part of implementing the main class is the constructor. It is unchanged from step-3 and step-4:
LaplaceProblem::LaplaceProblem () : fe (1), dof_handler (triangulation) {}
The following is the first function to be called in run()
. It sets up the mesh and degrees of freedom.
We start by creating the usual square mesh and changing the boundary indicator on the parts of the boundary where we have Dirichlet boundary conditions (top and bottom, i.e. faces two and three of the reference cell as defined by GeometryInfo), so that we can distinguish between the parts of the boundary where periodic and where Dirichlet boundary conditions hold. We then refine the mesh a fixed number of times, with child faces inheriting the boundary indicators previously set on the coarse mesh from their parents.
void LaplaceProblem::make_grid_and_dofs ()
{
GridGenerator::hyper_cube (triangulation);
triangulation.begin_active ()->face (2)->set_boundary_indicator (1);
triangulation.begin_active ()->face (3)->set_boundary_indicator (1);
triangulation.refine_global (5);
The next step is to distribute the degrees of freedom and produce a little bit of graphical output:
dof_handler.distribute_dofs (fe); std::cout << "Number of active cells: " << triangulation.n_active_cells () << std::endl << "Degrees of freedom: " << dof_handler.n_dofs () << std::endl;
Now it is the time for the constraints that come from the periodicity constraints. We do this in the following, separate function, after clearing any possible prior content from the constraints object:
constraints.clear (); make_periodicity_constraints ();
We also incorporate the homogeneous Dirichlet boundary conditions on the upper and lower parts of the boundary (i.e. the ones with boundary indicator 1) and close the ConstraintMatrix
object:
VectorTools::interpolate_boundary_values (dof_handler, 1, ZeroFunction<2> (), constraints); constraints.close ();
Then we create the sparsity pattern and the system matrix and initialize the solution and right-hand side vectors. This is again as in step-3 or step-6, for example:
CompressedSparsityPattern c_sparsity_pattern (dof_handler.n_dofs(),
dof_handler.n_dofs());
DoFTools::make_sparsity_pattern (dof_handler,
c_sparsity_pattern,
constraints,
false);
c_sparsity_pattern.compress ();
sparsity_pattern.copy_from (c_sparsity_pattern);
system_matrix.reinit (sparsity_pattern);
system_rhs.reinit (dof_handler.n_dofs());
solution.reinit (dof_handler.n_dofs());
}
This is the function that provides the new material of this tutorial program. The general outline of the algorithm is as follows: we first loop over all the degrees of freedom on the right boundary and record their -locations in a map together with their global indices. Then we go along the left boundary, find matching
-locations for each degree of freedom, and then add constraints that identify these matched degrees of freedom.
In this function, we make use of the fact that we have a scalar element (i.e. the only valid vector component that can be passed to DoFAccessor::vertex_dof_index is zero) and that we have a element for which all degrees of freedom live in the vertices of the cell. Furthermore, we have assumed that we are in 2d and that meshes were not refined adaptively — the latter assumption would imply that there may be vertices that aren't matched one-to-one and for which we won't be able to compute constraints this easily. We will discuss in the "outlook" part of the results section below other strategies to write the current function that can work in cases like this as well.
void LaplaceProblem::make_periodicity_constraints ()
{
To start with the actual implementation, we loop over all active cells and check whether the cell is located at the right boundary (i.e. face 1 — the one at the right end of the cell — is at the boundary). If that is so, then we use that for the currently used finite element, each degree of freedom of the face is located on one vertex, and store their -coordinate along with the global number of this degree of freedom in the following map:
std::map<unsigned int, double> dof_locations; for (DoFHandler<2>::active_cell_iterator cell = dof_handler.begin_active (); cell != dof_handler.end (); ++cell) if (cell->at_boundary () && cell->face(1)->at_boundary ()) { dof_locations[cell->face(1)->vertex_dof_index(0, 0)] = cell->face(1)->vertex(0)[1]; dof_locations[cell->face(1)->vertex_dof_index(1, 0)] = cell->face(1)->vertex(1)[1]; }
Note that in the above block, we add vertices zero and one of the affected face to the map. This means that we will add each vertex twice, once from each of the two adjacent cells (unless the vertex is a corner of the domain). Since the coordinates of the vertex are the same both times of course, there is no harm: we replace one value in the map with itself the second time we visit an entry.
The same will be true below where we add the same constraint twice to the ConstraintMatrix — again, we will overwrite the constraint with itself, and no harm is done.
Now we have to find the corresponding degrees of freedom on the left part of the boundary. Therefore we loop over all cells again and choose the ones where face 0 is at the boundary:
for (DoFHandler<2>::active_cell_iterator cell = dof_handler.begin_active (); cell != dof_handler.end (); ++cell) if (cell->at_boundary () && cell->face (0)->at_boundary ()) {
Every degree of freedom on this face needs to have a corresponding one on the right side of the face, and our goal is to add a constraint for the one on the left in terms of the one on the right. To this end we first add a new line to the constraint matrix for this one degree of freedom. Then we identify it with the corresponding degree of freedom on the right part of the boundary by constraining the degree of freedom on the left with the one on the right times a weight of 1.0.
Consequently, we loop over the two vertices of each face we find and then loop over all the -locations we've previously recorded to find which degree of freedom on the right boundary corresponds to the one we currently look at. Note that we have entered these into a map, and when looping over the iterators
p
of this map, p->first
corresponds to the "key" of an entry (the global number of the degree of freedom), whereas p->second
is the "value" (the -location we have entered above).
We are quite sure here that we should be finding such a corresponding degree of freedom. However, sometimes stuff happens and so the bottom of the block contains an assertion that our assumption was indeed correct and that a vertex was found.
for (unsigned int face_vertex = 0; face_vertex<2; ++face_vertex) { constraints.add_line (cell->face(0)->vertex_dof_index (face_vertex, 0)); std::map<unsigned int, double>::const_iterator p = dof_locations.begin(); for (; p != dof_locations.end(); ++p) if (std::fabs(p->second - cell->face(0)->vertex(face_vertex)[1]) < 1e-8) { constraints.add_entry (cell->face(0)->vertex_dof_index (face_vertex, 0), p->first, 1.0); break; } Assert (p != dof_locations.end(), ExcMessage ("No corresponding degree of freedom was found!")); } } }
Assembling the system matrix and the right-hand side vector is done as in other tutorials before.
The only difference here is that we don't copy elements from local contributions into the global matrix and later fix up constrained degrees of freedom, but that we let the ConstraintMatrix do this job in one swoop for us using the ConstraintMatrix::distribute_local_to_global function(). This was previously already demonstrated in step-16, step-22, for example, along with a discussion in the introduction of step-27.
void LaplaceProblem::assemble_system () { QGauss<2> quadrature_formula(2); FEValues<2> fe_values (fe, quadrature_formula, update_values | update_gradients | update_quadrature_points | update_JxW_values); const unsigned int dofs_per_cell = fe.dofs_per_cell; const unsigned int n_q_points = quadrature_formula.size(); FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell); Vector<double> cell_rhs (dofs_per_cell); std::vector<unsigned int> local_dof_indices (dofs_per_cell); const RightHandSide right_hand_side; DoFHandler<2>::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); for (; cell!=endc; ++cell) { fe_values.reinit (cell); cell_matrix = 0; cell_rhs = 0; for (unsigned int q_point=0; q_point<n_q_points; ++q_point) for (unsigned int i=0; i<dofs_per_cell; ++i) { for (unsigned int j=0; j<dofs_per_cell; ++j) cell_matrix(i,j) += (fe_values.shape_grad (i, q_point) * fe_values.shape_grad (j, q_point) * fe_values.JxW (q_point)); cell_rhs(i) += (fe_values.shape_value (i, q_point) * right_hand_side.value (fe_values.quadrature_point (q_point)) * fe_values.JxW (q_point)); } cell->get_dof_indices (local_dof_indices); constraints.distribute_local_to_global (cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs); } }
To solve the linear system of equations we use the CG solver with an SSOR-preconditioner. This is, again, copied almost verbatim from step-4, with the exception of the preconditioner. As in step-6, we need to make sure that constrained degrees of freedom get their correct values after solving by calling the ConstraintMatrix::distribute function:
void LaplaceProblem::solve ()
{
SolverControl solver_control (dof_handler.n_dofs (), 1e-12);
PreconditionSSOR<SparseMatrix<double> > precondition;
precondition.initialize (system_matrix);
SolverCG<> cg (solver_control);
cg.solve (system_matrix, solution, system_rhs, precondition);
constraints.distribute (solution);
}
This is another function copied from previous tutorial programs. It generates graphical output in VTK format:
void LaplaceProblem::output_results () { DataOut<2> data_out; data_out.attach_dof_handler (dof_handler); data_out.add_data_vector (solution, "u"); data_out.build_patches (); std::ofstream output ("solution.vtk"); data_out.write_vtk (output); }
And another function copied from previous programs:
void LaplaceProblem::run ()
{
make_grid_and_dofs();
assemble_system ();
solve ();
output_results ();
}
main
functionAnd at the end we have the main function as usual, this time copied from step-6:
int main () { try { deallog.depth_console (0); LaplaceProblem laplace_problem; laplace_problem.run (); } catch (std::exception &exc) { std::cerr << std::endl << std::endl << "----------------------------------------------------" << std::endl; std::cerr << "Exception on processing: " << std::endl << exc.what() << std::endl << "Aborting!" << std::endl << "----------------------------------------------------" << std::endl; return 1; } catch (...) { std::cerr << std::endl << std::endl << "----------------------------------------------------" << std::endl; std::cerr << "Unknown exception!" << std::endl << "Aborting!" << std::endl << "----------------------------------------------------" << std::endl; return 1; } return 0; }
The textual output the program generates is not very surprising. It just prints out the usual information on degrees of freedom and active cells, in much the same way as step-3 did:
Number of active cells: 1024 Degrees of freedom: 1089
The plot of the solution can be found in the figure below. We can see that the solution is constant zero on the upper and the lower part of the boundary as required by the homogeneous Dirichlet boundary conditions. On the left and right parts the values coincide with each other, just as we wanted:
Note also that the solution is clearly not left-right symmetric and so would not likely have been periodic had we prescribed, for example, homogeneous Neumann boundary condition. However, it is periodic thanks to the constraints imposed.
The function LaplaceProblem::make_periodicity_constraints
is relatively simple in that it just matches the location of degrees of freedom. This makes it flexible when the periodicity boundary conditions are posed not just on opposite faces of the unit rectangle but on separate parts of a possibly more complicated domain. Or, if one wanted to "twist" the boundary condition by prescribing, for example,
On the other hand, the function is somewhat limited by the assumption that the domain is two-dimensional and that we only use elements. The former assumption is easily lifted by looping over all four vertices of a face in 3d, but the latter is somewhat more complicated to lift because we have assumed that degrees of freedom are only located in vertices. In the following, therefore, let us describe a function that computes the same constraints but in a dimension-independent way and for any finite element one may want to consider.
The idea is to work recursively on pairs of faces. For example, let us start with the left and right face of the (single) coarse mesh cell. They need to match, but they are not active (i.e. they are further refined) and so there are no degrees of freedom associated with these faces. However, if the two current faces are periodic, then so are the zeroth children of the two as well as the respective first children, etc. We can then in turn work on each of these pairs of faces. If they are not active, we may recurse further into this refinement tree until we find a pair of active faces. In that case, we enter equivalences between matching degrees of freedom on the two active faces.
An implementation of this idea would look like follows (with the make_periodicity_constraint_recursively()
function — an implementation detail, not an external interface — put into an anonymous namespace):
namespace { template <int dim> void make_periodicity_constraints_recursively (const typename DoFHandler<dim>::face_iterator &face_1, const typename DoFHandler<dim>::face_iterator &face_2, ConstraintMatrix &constraints) { Assert (face_1->n_children() == face_2->n_children(), ExcNotImplemented()); if (face_1->has_children()) { for (unsigned int c=0; c<face_1->n_children(); ++c) make_periodicity_constraints_recursively<dim> (face_1->child(c), face_2->child(c), constraints); } else { const unsigned int dofs_per_face = face_1->get_dof_handler().get_fe().dofs_per_face; std::vector<unsigned int> local_dof_indices_1 (dofs_per_face); face_1->get_dof_indices (local_dof_indices_1); std::vector<unsigned int> local_dof_indices_2 (dofs_per_face); face_2->get_dof_indices (local_dof_indices_2); for (unsigned int i=0; i<dofs_per_face; ++i) { constraints.add_line (local_dof_indices_1[i]); constraints.add_entry (local_dof_indices_1[i], local_dof_indices_2[i], 1.0); } } } } void LaplaceProblem::make_periodicity_constraints () { make_periodicity_constraints_recursively<2> (dof_handler.begin(0)->face(0), dof_handler.begin(0)->face(1), constraints); }
The implementation of the recursive function should be mostly self explanatory given the discussion above. The LaplaceProblem::make_periodicity_constraints()
function simply calls the former with matching faces of the first (and only) coarse mesh cell on refinement level 0. Note that when calling the recursive function we have to explicitly specify the template argument since the compiler can not deduce it (the template argument is only used in a "non-deducible context").
This function is now dimension and finite element independent, but it still has the restriction that it assumes that the mesh is uniformly refined (or, in fact, only that matching periodic faces are refined equally). We check this at the beginning by asserting that both faces have the same number of children (that includes that neither have any children, i.e. that both are active). On the other hand, the function above can be extended to also allow this sort of thing. In that case, if we encounter the situation that only one cell is refined, we would have to recurse into its children and interpolate their degrees of freedom with respect to the degrees of freedom to the coarser matching face. This can use the same facilities the finite element classes already provide for computing constraints based on hanging nodes. We leave implementing this as an exercise, however.
(If you are looking at a locally installed deal.II version, then the program can be found at /build /buildd /deal.ii-6.3.1 /examples /step-45 /step-45.cc . Otherwise, this is only the path on some remote server.)
/ * Subversion Id: @ref step_45 "step-45".cc 21309 2010-06-24 03:02:55Z bangerth * / / * Author: Markus Buerg, University of Karlsruhe, 2010 * / / * Subversion Id: @ref step_45 "step-45".cc 21309 2010-06-24 03:02:55Z bangerth * / / * * / / * Copyright (C) 2010 by the deal.II authors * / / * * / / * This file is subject to QPL and may not be distributed * / / * without copyright and license information. Please refer * / / * to the file deal.II/doc/license.html for the text and * / / * further information on this license. * /
#include <base/function.h> #include <base/quadrature_lib.h> #include <lac/constraint_matrix.h> #include <lac/precondition.h> #include <lac/solver_cg.h> #include <lac/solver_control.h> #include <lac/sparse_matrix.h> #include <lac/sparsity_pattern.h> #include <lac/compressed_sparsity_pattern.h> #include <grid/grid_generator.h> #include <grid/tria.h> #include <dofs/dof_accessor.h> #include <dofs/dof_handler.h> #include <dofs/dof_tools.h> #include <fe/fe_q.h> #include <fe/fe_values.h> #include <numerics/data_out.h> #include <numerics/vectors.h> #include <fstream> using namespace dealii;
class LaplaceProblem { public: LaplaceProblem (); void run (); private: Triangulation<2> triangulation; FE_Q<2> fe; DoFHandler<2> dof_handler; ConstraintMatrix constraints; SparsityPattern sparsity_pattern; SparseMatrix<double> system_matrix; Vector<double> system_rhs; Vector<double> solution; void assemble_system (); void output_results (); void make_grid_and_dofs (); void make_periodicity_constraints (); void solve (); };
class RightHandSide: public Function<2> { public: RightHandSide (); virtual double value (const Point<2>& p, const unsigned int component = 0) const; }; RightHandSide::RightHandSide () : Function<2> () {} double RightHandSide::value (const Point<2>&p, const unsigned int) const { return (std::cos (2 * numbers::PI * p(0)) * std::exp (- 2 * p(0)) * std::cos (2 * numbers::PI * p(1)) * std::exp (- 2 * p(1))); }
LaplaceProblem::LaplaceProblem () : fe (1), dof_handler (triangulation) {}
void LaplaceProblem::make_grid_and_dofs () { GridGenerator::hyper_cube (triangulation); triangulation.begin_active ()->face (2)->set_boundary_indicator (1); triangulation.begin_active ()->face (3)->set_boundary_indicator (1); triangulation.refine_global (5); dof_handler.distribute_dofs (fe); std::cout << "Number of active cells: " << triangulation.n_active_cells () << std::endl << "Degrees of freedom: " << dof_handler.n_dofs () << std::endl; constraints.clear (); make_periodicity_constraints (); VectorTools::interpolate_boundary_values (dof_handler, 1, ZeroFunction<2> (), constraints); constraints.close (); CompressedSparsityPattern c_sparsity_pattern (dof_handler.n_dofs(), dof_handler.n_dofs()); DoFTools::make_sparsity_pattern (dof_handler, c_sparsity_pattern, constraints, false); c_sparsity_pattern.compress (); sparsity_pattern.copy_from (c_sparsity_pattern); system_matrix.reinit (sparsity_pattern); system_rhs.reinit (dof_handler.n_dofs()); solution.reinit (dof_handler.n_dofs()); }
void LaplaceProblem::make_periodicity_constraints () { std::map<unsigned int, double> dof_locations; for (DoFHandler<2>::active_cell_iterator cell = dof_handler.begin_active (); cell != dof_handler.end (); ++cell) if (cell->at_boundary () && cell->face(1)->at_boundary ()) { dof_locations[cell->face(1)->vertex_dof_index(0, 0)] = cell->face(1)->vertex(0)[1]; dof_locations[cell->face(1)->vertex_dof_index(1, 0)] = cell->face(1)->vertex(1)[1]; } for (DoFHandler<2>::active_cell_iterator cell = dof_handler.begin_active (); cell != dof_handler.end (); ++cell) if (cell->at_boundary () && cell->face (0)->at_boundary ()) { for (unsigned int face_vertex = 0; face_vertex<2; ++face_vertex) { constraints.add_line (cell->face(0)->vertex_dof_index (face_vertex, 0)); std::map<unsigned int, double>::const_iterator p = dof_locations.begin(); for (; p != dof_locations.end(); ++p) if (std::fabs(p->second - cell->face(0)->vertex(face_vertex)[1]) < 1e-8) { constraints.add_entry (cell->face(0)->vertex_dof_index (face_vertex, 0), p->first, 1.0); break; } Assert (p != dof_locations.end(), ExcMessage ("No corresponding degree of freedom was found!")); } } }
void LaplaceProblem::assemble_system () { QGauss<2> quadrature_formula(2); FEValues<2> fe_values (fe, quadrature_formula, update_values | update_gradients | update_quadrature_points | update_JxW_values); const unsigned int dofs_per_cell = fe.dofs_per_cell; const unsigned int n_q_points = quadrature_formula.size(); FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell); Vector<double> cell_rhs (dofs_per_cell); std::vector<unsigned int> local_dof_indices (dofs_per_cell); const RightHandSide right_hand_side; DoFHandler<2>::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); for (; cell!=endc; ++cell) { fe_values.reinit (cell); cell_matrix = 0; cell_rhs = 0; for (unsigned int q_point=0; q_point<n_q_points; ++q_point) for (unsigned int i=0; i<dofs_per_cell; ++i) { for (unsigned int j=0; j<dofs_per_cell; ++j) cell_matrix(i,j) += (fe_values.shape_grad (i, q_point) * fe_values.shape_grad (j, q_point) * fe_values.JxW (q_point)); cell_rhs(i) += (fe_values.shape_value (i, q_point) * right_hand_side.value (fe_values.quadrature_point (q_point)) * fe_values.JxW (q_point)); } cell->get_dof_indices (local_dof_indices); constraints.distribute_local_to_global (cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs); } }
void LaplaceProblem::solve ()
{
SolverControl solver_control (dof_handler.n_dofs (), 1e-12);
PreconditionSSOR<SparseMatrix<double> > precondition;
precondition.initialize (system_matrix);
SolverCG<> cg (solver_control);
cg.solve (system_matrix, solution, system_rhs, precondition);
constraints.distribute (solution);
}
void LaplaceProblem::output_results () { DataOut<2> data_out; data_out.attach_dof_handler (dof_handler); data_out.add_data_vector (solution, "u"); data_out.build_patches (); std::ofstream output ("solution.vtk"); data_out.write_vtk (output); }
void LaplaceProblem::run ()
{
make_grid_and_dofs();
assemble_system ();
solve ();
output_results ();
}
int main () { try { deallog.depth_console (0); LaplaceProblem laplace_problem; laplace_problem.run (); } catch (std::exception &exc) { std::cerr << std::endl << std::endl << "----------------------------------------------------" << std::endl; std::cerr << "Exception on processing: " << std::endl << exc.what() << std::endl << "Aborting!" << std::endl << "----------------------------------------------------" << std::endl; return 1; } catch (...) { std::cerr << std::endl << std::endl << "----------------------------------------------------" << std::endl; std::cerr << "Unknown exception!" << std::endl << "Aborting!" << std::endl << "----------------------------------------------------" << std::endl; return 1; } return 0; }