void
PointValueEvaluation<dim>::
operator () (const DoFHandler<dim> &dof_handler,
- const Vector<double> &solution) const
+ const Vector<double> &solution) const
{
- // First allocate a variable that
- // will hold the point
- // value. Initialize it with a
- // value that is clearly bogus,
- // so that if we fail to set it
- // to a reasonable value, we will
- // note at once. This may not be
- // necessary in a function as
- // small as this one, since we
- // can easily see all possible
- // paths of execution here, but
- // it proved to be helpful for
- // more complex cases, and so we
- // employ this strategy here as
- // well.
+ // First allocate a variable that will hold the point value. Initialize
+ // it with a value that is clearly bogus, so that if we fail to set it
+ // to a reasonable value, we will note at once. This may not be
+ // necessary in a function as small as this one, since we can easily see
+ // all possible paths of execution here, but it proved to be helpful for
+ // more complex cases, and so we employ this strategy here as well.
double point_value = 1e20;
- // Then loop over all cells and
- // all their vertices, and check
- // whether a vertex matches the
- // evaluation point. If this is
- // the case, then extract the
- // point value, set a flag that
- // we have found the point of
+ // Then loop over all cells and all their vertices, and check whether a
+ // vertex matches the evaluation point. If this is the case, then
+ // extract the point value, set a flag that we have found the point of
// interest, and exit the loop.
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
void
PointXDerivativeEvaluation<dim>::
operator () (const DoFHandler<dim> &dof_handler,
- const Vector<double> &solution) const
+ const Vector<double> &solution) const
{
- // This time initialize the
- // return value with something
- // useful, since we will have to
- // add up a number of
- // contributions and take the
- // mean value afterwards...
+ // This time initialize the return value with something useful, since we
+ // will have to add up a number of contributions and take the mean value
+ // afterwards...
double point_derivative = 0;
- // ...then have some objects of
- // which the meaning wil become
- // clear below...
+ // ...then have some objects of which the meaning wil become clear
+ // below...
QTrapez<dim> vertex_quadrature;
FEValues<dim> fe_values (dof_handler.get_fe(),
vertex_quadrature,
};
- // And then we also have to define
- // these respective functions, of
- // course. Given our discussion in
- // the introduction of how the
- // solution should look like, the
- // following computations should be
- // straightforward:
+ // And then we also have to define these respective functions, of
+ // course. Given our discussion in the introduction of how the solution
+ // should look like, the following computations should be straightforward:
template <int dim>
- double RightHandSide<dim>::value (const Point<dim> & /*p*/,
+ double RightHandSide<dim>::value (const Point<dim> &/*p*/,
const unsigned int /*component*/) const
{
return 0;
// @sect5{Private member variables}
- // Next, we have a few member
- // variables. In particular,
- // these are (i) a reference to
- // the parameter object (owned by
- // the main function of this
- // program, and passed to the
- // constructor of this class),
- // (ii) an object describing the
- // material parameters for the
- // number of energy groups
- // requested in the input file,
- // and (iii) the finite element
- // to be used by all energy
- // groups:
+ // Next, we have a few member variables. In particular, these are (i) a
+ // reference to the parameter object (owned by the main function of this
+ // program, and passed to the constructor of this class), (ii) an object
+ // describing the material parameters for the number of energy groups
+ // requested in the input file, and (iii) the finite element to be used by
+ // all energy groups:
- const Parameters ¶meters;
+ const Parameters ¶meters;
const MaterialData material_data;
FE_Q<dim> fe;
- // The constructor takes the
- // ParameterHandler object and stores
- // it in a reference. It also
- // initializes the DoF-Handler and
- // the finite element system, which
- // consists of two copies of the
- // scalar Q1 field, one for $v$ and
- // one for $w$:
+ // The constructor takes the ParameterHandler object and stores it in a
+ // reference. It also initializes the DoF-Handler and the finite element
+ // system, which consists of two copies of the scalar Q1 field, one for $v$
+ // and one for $w$:
template <int dim>
- UltrasoundProblem<dim>::UltrasoundProblem (ParameterHandler ¶m)
+ UltrasoundProblem<dim>::UltrasoundProblem (ParameterHandler ¶m)
:
prm(param),
dof_handler(triangulation),
EulerEquations<dim>::Postprocessor::
compute_derived_quantities_vector (const std::vector<Vector<double> > &uh,
const std::vector<std::vector<Tensor<1,dim> > > &duh,
- const std::vector<std::vector<Tensor<2,dim> > > & /*dduh*/,
- const std::vector<Point<dim> > & /*normals*/,
- const std::vector<Point<dim> > & /*evaluation_points*/,
+ const std::vector<std::vector<Tensor<2,dim> > > &/*dduh*/,
+ const std::vector<Point<dim> > &/*normals*/,
+ const std::vector<Point<dim> > &/*evaluation_points*/,
std::vector<Vector<double> > &computed_quantities) const
{
- // At the beginning of the function, let us
- // make sure that all variables have the
- // correct sizes, so that we can access
- // individual vector elements without
- // having to wonder whether we might read
- // or write invalid elements; we also check
- // that the <code>duh</code> vector only
- // contains data if we really need it (the
- // system knows about this because we say
- // so in the
- // <code>get_needed_update_flags()</code>
- // function below). For the inner vectors,
- // we check that at least the first element
- // of the outer vector has the correct
- // inner size:
+ // At the beginning of the function, let us make sure that all variables
+ // have the correct sizes, so that we can access individual vector
+ // elements without having to wonder whether we might read or write
+ // invalid elements; we also check that the <code>duh</code> vector only
+ // contains data if we really need it (the system knows about this because
+ // we say so in the <code>get_needed_update_flags()</code> function
+ // below). For the inner vectors, we check that at least the first element
+ // of the outer vector has the correct inner size:
const unsigned int n_quadrature_points = uh.size();
if (do_schlieren_plot == true)
// @sect3{Include files}
//
- // Most of the include files we need for this
- // program have already been discussed in
- // previous programs. In particular, all of
- // the following should already be familiar
- // friends:
+ // Most of the include files we need for this program have already been
+ // discussed in previous programs. In particular, all of the following should
+ // already be familiar friends:
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
+
+#include <deal.II/lac/abstract_linear_algebra.h>
+
+#define USE_PETSC_LA
+
+namespace LA
+{
+#ifdef USE_PETSC_LA
+ using namespace dealii::LinearAlgebraPETSc;
+#else
+ using namespace dealii::LinearAlgebraTrilinos;
+#endif
+}
+
+
#include <deal.II/lac/vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/solver_cg.h>
locally_relevant_dofs);
locally_relevant_solution = 0;
system_rhs.reinit (mpi_communicator,
- dof_handler.n_dofs(),
- dof_handler.n_locally_owned_dofs());
+ locally_owned_dofs);
+
system_rhs = 0;
- // The next step is to compute hanging node
- // and boundary value constraints, which we
- // combine into a single object storing all
+ // The next step is to compute hanging node and boundary value
+ // constraints, which we combine into a single object storing all
// constraints.
//
- // As with all other things in %parallel,
- // the mantra must be that no processor can
- // store all information about the entire
- // universe. As a consequence, we need to
- // tell the constraints object for which
- // degrees of freedom it can store
- // constraints and for which it may not
- // expect any information to store. In our
- // case, as explained in the @ref
- // distributed module, the degrees of
- // freedom we need to care about on each
- // processor are the locally relevant ones,
- // so we pass this to the
- // ConstraintMatrix::reinit function. As a
- // side note, if you forget to pass this
- // argument, the ConstraintMatrix class
- // will allocate an array with length equal
- // to the largest DoF index it has seen so
- // far. For processors with high MPI
- // process number, this may be very large
- // -- maybe on the order of billions. The
- // program would then allocate more memory
- // than for likely all other operations
- // combined for this single array.
+ // As with all other things in %parallel, the mantra must be that no
+ // processor can store all information about the entire universe. As a
+ // consequence, we need to tell the constraints object for which degrees
+ // of freedom it can store constraints and for which it may not expect any
+ // information to store. In our case, as explained in the @ref distributed
+ // module, the degrees of freedom we need to care about on each processor
+ // are the locally relevant ones, so we pass this to the
+ // ConstraintMatrix::reinit function. As a side note, if you forget to
+ // pass this argument, the ConstraintMatrix class will allocate an array
+ // with length equal to the largest DoF index it has seen so far. For
+ // processors with high MPI process number, this may be very large --
+ // maybe on the order of billions. The program would then allocate more
+ // memory than for likely all other operations combined for this single
+ // array.
constraints.clear ();
constraints.reinit (locally_relevant_dofs);
DoFTools::make_hanging_node_constraints (dof_handler, constraints);
SolverControl solver_control (dof_handler.n_dofs(), 1e-12);
+#ifdef USE_PETSC_LA
+
PETScWrappers::SolverCG solver(solver_control, mpi_communicator);
- // Ask for a symmetric preconditioner by
- // setting the first parameter in
+ // Ask for a symmetric preconditioner by setting the first parameter in
// AdditionalData to true.
- PETScWrappers::PreconditionBoomerAMG
+ LA::MPI::PreconditionAMG
preconditioner(system_matrix,
PETScWrappers::PreconditionBoomerAMG::AdditionalData(true));
fe_values.JxW(q_point));
}
+ // Finally, transfer the contributions from @p cell_matrix and
+ // @p cell_rhs into the global objects.
cell->get_dof_indices (local_dof_indices);
- // transfer the contributions from @p cell_matrix and @cell_rhs into the global objects.
- constraints.distribute_local_to_global(cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
+ constraints.distribute_local_to_global(cell_matrix,
+ cell_rhs,
+ local_dof_indices,
+ system_matrix, system_rhs);
}
- // Now we are done assembling the linear
- // system. The constrained nodes are still
- // in the linear system (there is a one on
- // the diagonal of the matrix and all other
- // entries for this line are set to zero)
- // but the computed values are invalid. We
- // compute the correct values for these
- // nodes at the end of the
- // <code>solve</code> function.
+
+#ifdef USE_PETSC_LA
+ system_matrix.compress();
+ system_rhs.compress();
+#endif
+
+ // Now we are done assembling the linear system. The constrained nodes are
+ // still in the linear system (there is a one on the diagonal of the matrix
+ // and all other entries for this line are set to zero) but the computed
+ // values are invalid. We compute the correct values for these nodes at the
+ // end of the <code>solve</code> function.
}
// @sect3{The <code>main</code> function}
- // The main function is unaltered in
- // its functionality from the
- // previous example, but we have
- // taken a step of additional
- // caution. Sometimes, something goes
- // wrong (such as insufficient disk
- // space upon writing an output file,
- // not enough memory when trying to
- // allocate a vector or a matrix, or
- // if we can't read from or write to
- // a file for whatever reason), and
- // in these cases the library will
- // throw exceptions. Since these are
- // run-time problems, not programming
- // errors that can be fixed once and
- // for all, this kind of exceptions
- // is not switched off in optimized
- // mode, in contrast to the
- // <code>Assert</code> macro which we have
- // used to test against programming
- // errors. If uncaught, these
- // exceptions propagate the call tree
- // up to the <code>main</code> function, and
- // if they are not caught there
- // either, the program is aborted. In
- // many cases, like if there is not
- // enough memory or disk space, we
- // can't do anything but we can at
- // least print some text trying to
- // explain the reason why the program
- // failed. A way to do so is shown in
- // the following. It is certainly
- // useful to write any larger program
- // in this way, and you can do so by
- // more or less copying this function
- // except for the <code>try</code> block that
- // actually encodes the functionality
- // particular to the present
- // application.
+ // The main function is unaltered in its functionality from the previous
+ // example, but we have taken a step of additional caution. Sometimes,
+ // something goes wrong (such as insufficient disk space upon writing an
+ // output file, not enough memory when trying to allocate a vector or a
+ // matrix, or if we can't read from or write to a file for whatever reason),
+ // and in these cases the library will throw exceptions. Since these are
+ // run-time problems, not programming errors that can be fixed once and for
+ // all, this kind of exceptions is not switched off in optimized mode, in
+ // contrast to the <code>Assert</code> macro which we have used to test
+ // against programming errors. If uncaught, these exceptions propagate the
+ // call tree up to the <code>main</code> function, and if they are not caught
+ // there either, the program is aborted. In many cases, like if there is not
+ // enough memory or disk space, we can't do anything but we can at least print
+ // some text trying to explain the reason why the program failed. A way to do
+ // so is shown in the following. It is certainly useful to write any larger
+ // program in this way, and you can do so by more or less copying this
+ // function except for the <code>try</code> block that actually encodes the
+ // functionality particular to the present application.
-int main ()
+int main (int argc, char *argv[])
{
- // The general idea behind the
- // layout of this function is as
- // follows: let's try to run the
- // program as we did before...
+ // The general idea behind the layout of this function is as follows: let's
+ // try to run the program as we did before...
try
{
deallog.depth_console (0);
+ Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv);
- Step6<2> laplace_problem_2d;
- laplace_problem_2d.run ();
+ {
+ Step6<2> laplace_problem_2d;
+ laplace_problem_2d.run ();
+ }
}
- // ...and if this should fail, try
- // to gather as much information as
- // possible. Specifically, if the
- // exception that was thrown is an
- // object of a class that is
- // derived from the C++ standard
- // class <code>exception</code>, then we can
- // use the <code>what</code> member function
- // to get a string which describes
- // the reason why the exception was
+ // ...and if this should fail, try to gather as much information as
+ // possible. Specifically, if the exception that was thrown is an object of
+ // a class that is derived from the C++ standard class
+ // <code>exception</code>, then we can use the <code>what</code> member
+ // function to get a string which describes the reason why the exception was
// thrown.
//
- // The deal.II exception classes
- // are all derived from the
- // standard class, and in
- // particular, the <code>exc.what()</code>
- // function will return
- // approximately the same string as
- // would be generated if the
- // exception was thrown using the
- // <code>Assert</code> macro. You have seen
- // the output of such an exception
- // in the previous example, and you
- // then know that it contains the
- // file and line number of where
- // the exception occured, and some
- // other information. This is also
- // what the following statements
- // would print.
+ // The deal.II exception classes are all derived from the standard class,
+ // and in particular, the <code>exc.what()</code> function will return
+ // approximately the same string as would be generated if the exception was
+ // thrown using the <code>Assert</code> macro. You have seen the output of
+ // such an exception in the previous example, and you then know that it
+ // contains the file and line number of where the exception occured, and
+ // some other information. This is also what the following statements would
+ // print.
//
- // Apart from this, there isn't
- // much that we can do except
- // exiting the program with an
- // error code (this is what the
- // <code>return 1;</code> does):
+ // Apart from this, there isn't much that we can do except exiting the
+ // program with an error code (this is what the <code>return 1;</code>
+ // does):
catch (std::exception &exc)
{
std::cerr << std::endl << std::endl
template <int dim>
void
GradientEstimation::estimate (const DoFHandler<dim> &dof_handler,
- const Vector<double> &solution,
+ const Vector<double> &solution,
Vector<float> &error_per_cell)
{
- // Before starting with the work,
- // we check that the vector into
- // which the results are written,
- // has the right size. It is a
- // common error that such
- // parameters have the wrong size,
- // but the resulting damage by not
- // catching these errors are very
- // subtle as they are usually
- // corruption of data somewhere in
- // memory. Often, the problems
- // emerging from this are not
- // reproducible, and we found that
- // it is well worth the effort to
+ // Before starting with the work, we check that the vector into which the
+ // results are written, has the right size. It is a common error that such
+ // parameters have the wrong size, but the resulting damage by not
+ // catching these errors are very subtle as they are usually corruption of
+ // data somewhere in memory. Often, the problems emerging from this are
+ // not reproducible, and we found that it is well worth the effort to
// check for such things.
Assert (error_per_cell.size() == dof_handler.get_tria().n_active_cells(),
ExcInvalidVectorLength (error_per_cell.size(),
/**
* Name of the exception and call sequence.
*/
- const char *exc;
+ const char *exc;
/**
- * A backtrace to the position
- * where the problem happened, if
- * the system supports this.
+ * A backtrace to the position where the problem happened, if the
+ * system supports this.
*/
char **stacktrace;
projects / dealii-svn.git / commitdiff