// of value 1 over the whole
// computational domain, i.e. by
// computing the areas $\int_K 1
- // dx=\int_{\hat K} 1 J(\hat x) d\hat
- // x \approx \sum J(\hat x)w(\hat x)$
- // of all active cells of the
- // triangulation and summing up these
- // contributions to gain the area of
- // the overall domain. The integrals
- // on each cell are approximated by
- // numerical quadrature, hence the
- // only additional ingredient we need
- // is to set up a FEValues object
- // that provides the corresponding
- // `JxW' values of each cell. We note
- // that here we won't use the
+ // dx=\int_{\hat K} 1 det J(\hat x)
+ // d\hat x \approx \sum det J(\hat
+ // xi)w(\hat xi)$, where the sum
+ // extends over all quadrature points
+ // on all active cells in the
+ // triangulation, with $w(xi)$ being
+ // the weight of quadrature point
+ // $xi$. The integrals on each cell
+ // are approximated by numerical
+ // quadrature, hence the only
+ // additional ingredient we need is
+ // to set up a FEValues object that
+ // provides the corresponding `JxW'
+ // values of each cell. (Note that
+ // `JxW' is meant to abbreviate
+ // ``Jacobian determinant times
+ // weight''; since in numerical
+ // quadrature the two factors always
+ // occur at the same places, we only
+ // offer the combined quantity,
+ // rather than two separate ones.) We
+ // note that here we won't use the
// FEValues object in its original
- // purpose that is computing the
- // values of basis functions of a
- // specific finite element. But here
- // we use it only to gain the `JxW'
- // at the quadrature points,
- // irrespective of the (dummy) finite
- // element we will give to the
- // constructor of the FEValues
- // object.
+ // purpose, i.e. for the computation
+ // of values of basis functions of a
+ // specific finite element at certain
+ // quadrature points. Rather, we use
+ // it only to gain the `JxW' at the
+ // quadrature points, irrespective of
+ // the (dummy) finite element we will
+ // give to the constructor of the
+ // FEValues object. The actual finite
+ // element given to the FEValues
+ // object is not used at all, so we
+ // could give any.
template <int dim>
void compute_pi_by_area ()
{
// For the numerical quadrature on
// all cells we employ a quadrature
// rule of sufficiently high
- // degree. We choose QGauss4 that is
- // of order 8, to be sure that the
- // errors due to numerical
+ // degree. We choose QGauss4 that
+ // is of order 8, to be sure that
+ // the errors due to numerical
// quadrature are of higher order
// than the order (maximal 6) that
// will occur due to the order of
// the approximation of the
// boundary, i.e. the order of the
- // mappings employed.
+ // mappings employed. Note that the
+ // integrand, the Jacobian
+ // determinant, is not a polynomial
+ // function (rather, it is a
+ // rational one), so we do not use
+ // Gauss quadrature in order to get
+ // the exact value of the integral
+ // as done often in finite element
+ // computations, but could as well
+ // have used any quadrature formula
+ // of like order instead.
const QGauss4<dim> quadrature;
// Now start by looping over
- // degrees=1..4
+ // polynomial mapping degrees=1..4:
for (unsigned int degree=1; degree<5; ++degree)
{
std::cout << "Degree = " << degree << std::endl;
- // Then we generate the
+ // First generate the
// triangulation, the Boundary
// and the Mapping object as
// already seen.
// We now create a dummy finite
// element. Here we could
- // choose a finite element no
- // matter which, as we are only
- // interested in the `JxW'
- // values provided by the
- // FEValues object below.
+ // choose any finite element,
+ // as we are only interested in
+ // the `JxW' values provided by
+ // the FEValues object
+ // below. Nevertheless, we have
+ // to provide a finite element
+ // since in this example we
+ // abuse the FEValues class a
+ // little in that we only ask
+ // it to provide us with the
+ // weights of certain
+ // quadrature points, in
+ // contrast to the usual
+ // purpose (and name) of the
+ // FEValues class which is to
+ // provide the values of finite
+ // elements at these points.
const FE_Q<dim> dummy_fe (1);
- // Then we create a DofHandler
- // object. This object will
- // provide us with
+ // Likewise, we need to create
+ // a DofHandler object. We do
+ // not actually use it, but it
+ // will provide us with
// `active_cell_iterators' that
// are needed to reinit the
// FEValues object on each cell
// the dummy finite element and
// the quadrature object to the
// constructor, together with
- // the UpdateFlag asking for
+ // the UpdateFlags asking for
// the `JxW' values at the
- // quadrature points only.
+ // quadrature points only. This
+ // tells the FEValues object
+ // that it needs not compute
+ // other quantities upon
+ // calling the ``reinit''
+ // function, thus saving
+ // computation time.
+ //
+ // The most important
+ // difference in the
+ // construction of the FEValues
+ // object compared to previous
+ // example programs is that we
+ // pass a mapping object as
+ // first argument, which is to
+ // be used in the computation
+ // of the mapping from unit to
+ // real cell. In previous
+ // examples, this argument was
+ // omitted, resulting in the
+ // implicit use of an object of
+ // type MappingQ1.
FEValues<dim> fe_values (mapping, dummy_fe, quadrature, update_JxW_values);
// We employ an object of the
// ConvergenceTable class to
// store all important data
- // like the approximative
- // values for pi and the error
- // wrt. the true value of
- // pi. We will use functions
- // provided by the
+ // like the approximated values
+ // for pi and the error with
+ // respect to the true value of
+ // pi. We will also use
+ // functions provided by the
// ConvergenceTable class to
// compute convergence rates of
// the approximations to pi.
// automatically creates a
// table column with
// superscription `cells',
- // for the case this column
- // was not created before.
+ // in case this column was
+ // not created before.
table.add_value("cells", triangulation.n_active_cells());
// Then we distribute the
- // degrees of freedoms for
+ // degrees of freedom for
// the dummy finite
// element. Strictly
// speaking we do not need
long double area = 0;
// Now we loop over all
- // cells, reinit the
+ // cells, reinitialize the
// FEValues object for each
- // cell, add all `JxW'
- // values to `area'
+ // cell, and add up all the
+ // `JxW' values for this
+ // cell to `area'...
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
area += fe_values.JxW (i);
};
- // and store the resulting
- // area values and the
- // errors in the table. We
- // need a static cast to
- // double as there is no
- // add_value(string, long
- // double) function
+ // ...and store the
+ // resulting area values
+ // and the errors in the
+ // table. We need a static
+ // cast to double as there
+ // is no add_value(string,
+ // long double) function
// implemented.
table.add_value("eval.pi", static_cast<double> (area));
table.add_value("error", fabs(area-pi));
// `evaluate_all_convergence_rates'
table.omit_column_from_convergence_rate_evaluation("cells");
table.omit_column_from_convergence_rate_evaluation("eval.pi");
- table.evaluate_all_convergence_rates(
- ConvergenceTable::reduction_rate_log2);
+ table.evaluate_all_convergence_rates(ConvergenceTable::reduction_rate_log2);
// Finally we set the precision
- // and the scientific mode
+ // and scientific mode for
+ // output of some of the
+ // quantities...
table.set_precision("eval.pi", 16);
table.set_scientific("error", true);
- // and write the whole table to
- // cout.
+ // ...and write the whole table
+ // to std::cout.
table.write_text(std::cout);
std::cout << std::endl;
};
- // The following function also
+ // The following, second function also
// computes an approximation of pi
- // but this time via the diameter
+ // but this time via the perimeter
// 2*pi*radius of the domain instead
// of the area. This function is only
// a variation of the previous
DoFHandler<dim> dof_handler (triangulation);
- // Then we create a FEFaceValues
- // object instead of a FEValues
- // object as in the previous
- // function.
+ // Then we create a
+ // FEFaceValues object instead
+ // of a FEValues object as in
+ // the previous
+ // function. Again, we pass a
+ // mapping as first argument.
FEFaceValues<dim> fe_face_values (mapping, fe, quadrature, update_JxW_values);
ConvergenceTable table;
for (unsigned int i=0; i<fe_face_values.n_quadrature_points; ++i)
perimeter += fe_face_values.JxW (i);
};
- // We store the evaluated
- // values in the table
+ // Then store the evaluated
+ // values in the table...
table.add_value("eval.pi", static_cast<double> (perimeter/2.));
table.add_value("error", fabs(perimeter/2.-pi));
};
- // and we end this function as
- // we did in the previous
- // function.
+ // ...and end this function as
+ // we did in the previous one:
table.omit_column_from_convergence_rate_evaluation("cells");
table.omit_column_from_convergence_rate_evaluation("eval.pi");
- table.evaluate_all_convergence_rates(
- ConvergenceTable::reduction_rate_log2);
+ table.evaluate_all_convergence_rates(ConvergenceTable::reduction_rate_log2);
table.set_precision("eval.pi", 16);
table.set_scientific("error", true);
// The following main function just
- // calles the above functions in the
- // order of appearance.
+ // calls the above functions in the
+ // order of their appearance.
int main ()
{
std::cout.precision (16);
projects / dealii-svn.git / commitdiff