--- Begin Message ---
------------------------------------------------------------
revno: 5211
committer: Marie E. Rognes <meg@xxxxxxxxx>
branch nick: dolfin-error-control
timestamp: Fri 2010-11-05 13:06:31 +0100
message:
Add documentation (I think) and use the peculiar dolfin convention for
braces. Will continue with other adaptivity classes.
modified:
dolfin/adaptivity/ErrorControl.cpp
dolfin/adaptivity/ErrorControl.h
--
lp:~dolfin-core/dolfin/error-control
https://code.launchpad.net/~dolfin-core/dolfin/error-control
Your team DOLFIN Core Team is subscribed to branch lp:~dolfin-core/dolfin/error-control.
To unsubscribe from this branch go to https://code.launchpad.net/~dolfin-core/dolfin/error-control/+edit-subscription
=== modified file 'dolfin/adaptivity/ErrorControl.cpp'
--- dolfin/adaptivity/ErrorControl.cpp 2010-11-04 18:01:46 +0000
+++ dolfin/adaptivity/ErrorControl.cpp 2010-11-05 12:06:31 +0000
@@ -20,31 +20,13 @@
ErrorControl::ErrorControl(const Form& a, const Form& L,
const GoalFunctional& M)
: a(a), L(L), M(M)
-
{
- // The sub-classed constructor should initialize all forms AND
- // attach all relevant coefficients.
+ // Do nothing else
}
//-----------------------------------------------------------------------------
double ErrorControl::estimate_error(const Function& u,
- std::vector<const BoundaryCondition*> bcs) {
- /*
- w_h = E z_h - Pi_h E z_h,
-
- The error estimate eta_h is then:
-
- eta_h = r(w_h)
-
- FIXMEs:
-
- (1) Need to deal with possible sub-spaces for boundary condition
- spaces
-
- (2) Need to deal with casts more robustly. (Assume DirichletBC
- for now.)
-
- */
-
+ std::vector<const BoundaryCondition*> bcs)
+{
begin("Estimating error");
// Compute dual
@@ -57,40 +39,39 @@
// Attach dual and primal approximation (if not already attached) to
// residual
_residual->set_coefficient("__improved_dual", *_Ez_h);
- try {
+ try
+ {
_residual->set_coefficient("__discrete_primal_solution", u);
- } catch(...) {
+ } catch(...)
+ {
info("Cannot set coefficient named '__discrete_primal_solution'. Don't worry.");
}
// Assemble error estimate
const double error_estimate = assemble(*_residual);
-
end();
// Return estimate
return error_estimate;
-
}
//-----------------------------------------------------------------------------
void ErrorControl::compute_dual(Function& z,
- std::vector<const BoundaryCondition*> bcs) {
-
- // Create dual (homogenized) boundary conditions
+ std::vector<const BoundaryCondition*> bcs)
+{
+ // FIXME: Create dual (homogenized) boundary conditions
std::vector<const BoundaryCondition*> dual_bcs;
for (uint i = 0; i < bcs.size(); i++)
- dual_bcs.push_back(bcs[i]); // FIXME: push_back(bcs[i].homogenize());
+ dual_bcs.push_back(bcs[i]); // push_back(bcs[i].homogenize());
- // Create dual variational problem
+ // Create and solve dual variational problem
VariationalProblem dual(*_a_star, *_L_star, dual_bcs);
-
- // Solve dual problem
dual.solve(z);
}
//-----------------------------------------------------------------------------
void ErrorControl::compute_extrapolation(const Function& z,
- std::vector<const BoundaryCondition*> bcs){
-
+ std::vector<const BoundaryCondition*> bcs)
+{
+ // Extrapolate z
_Ez_h.reset(new Function(*_E));
_Ez_h->extrapolate(z);
@@ -98,18 +79,16 @@
Function u0(_Ez_h->function_space());
for (uint i = 0; i < bcs.size(); i++)
{
-
- // FIXME: Cast bc
+ // FIXME: Suboptimal cast.
DirichletBC bc(*dynamic_cast<const DirichletBC*>(bcs[i]));
+ // Extract SubSpace component
const FunctionSpace& V(bc.function_space());
-
- // Extract SubSpace component
const Array<uint>& component(V.component());
// If bcs[i].function_space is non subspace:
- if (component.size() == 0) {
-
+ if (component.size() == 0)
+ {
// Define constant 0.0 on this space
Function u0(V);
@@ -118,7 +97,6 @@
// Apply boundary condition to extrapolation
e_bc.apply(_Ez_h->vector());
-
continue;
}
@@ -136,14 +114,14 @@
}
}
//-----------------------------------------------------------------------------
-void ErrorControl::compute_indicators(Vector& indicators, const Function& u) {
-
+void ErrorControl::compute_indicators(Vector& indicators, const Function& u)
+{
begin("Computing indicators from error control");
- // R_T should live in _a_R_T's trial space while R_dT should live in
- // a special space
+ // Create Function for the strong cell residual (R_T)
Function R_T(*(_a_R_T->function_space(1)));
+ // Create SpecialFacetFunction for the strong facet residual (R_dT)
std::vector<Function *> f_e;
for (uint i=0; i <= R_T.geometric_dimension(); i++)
f_e.push_back(new Function(*(_a_R_dT->function_space(1))));
@@ -153,7 +131,8 @@
R_dT = new SpecialFacetFunction(f_e);
else if (f_e[0]->value_rank() == 1)
R_dT = new SpecialFacetFunction(f_e, f_e[0]->value_dimension(0));
- else {
+ else
+ {
R_dT = new SpecialFacetFunction(f_e, f_e[0]->value_dimension(0));
error("Not implemented for tensor-valued functions");
}
@@ -161,14 +140,13 @@
// Compute residual representation
residual_representation(R_T, *R_dT, u);
- // Check that extrapolation is computed (or compute it)
+ // FIXME: Check that extrapolation is computed (or compute it)
/// Interpolate dual extrapolation into primal test space
Function Pi_E_z_h(*(a.function_space(0)));
Pi_E_z_h.interpolate(*_Ez_h);
// Attach coefficients to error indicator form
- // NB: This depends on generated code order.
_eta_T->set_coefficient("__cell_residual", R_T);
_eta_T->set_coefficient("__facet_residual", *R_dT);
_eta_T->set_coefficient("__improved_dual", *_Ez_h);
@@ -177,7 +155,7 @@
// Assemble error indicator form
assemble(indicators, *_eta_T);
- // Take absolute value of indicators FIXME. Must be better way.
+ // Take absolute value of indicators: FIXME. Must be better way.
double abs;
for (uint i=0; i < indicators.size(); i++)
{
@@ -191,13 +169,12 @@
delete R_dT;
end();
-
}
//-----------------------------------------------------------------------------
void ErrorControl::residual_representation(Function& R_T,
SpecialFacetFunction& R_dT,
- const Function& u) {
-
+ const Function& u)
+{
begin("Computing residual representation");
// Compute cell residual
@@ -207,24 +184,24 @@
compute_facet_residual(R_dT, u, R_T);
end();
-
}
//-----------------------------------------------------------------------------
-void ErrorControl::compute_cell_residual(Function& R_T, const Function& u) {
-
- // Assuming that all coefficients except u_h are initialized prior
- // to this
+void ErrorControl::compute_cell_residual(Function& R_T, const Function& u)
+{
begin("Computing cell residual representation");
- // Attach u_h to _L_R_T
- try {
+ // Attach primal approximation to left-hand side form (residual) if
+ // necessary
+ try
+ {
_L_R_T->set_coefficient("__discrete_primal_solution", u);
- } catch(...) {
+ } catch(...)
+ {
info("Cannot set coefficient named 'the_discrete_solution'. Don't worry.");
}
- // Compute residual decomposition (FIXME: do this locally)
+ // Define variational problem for cell residual and solve
VariationalProblem cell_residual(*_a_R_T, *_L_R_T);
cell_residual.parameters["linear_solver"] = "iterative";
cell_residual.solve(R_T);
@@ -234,12 +211,10 @@
//-----------------------------------------------------------------------------
void ErrorControl::compute_facet_residual(SpecialFacetFunction& R_dT,
const Function& u,
- const Function& R_T) {
-
+ const Function& R_T)
+{
begin("Computing facet residual representation");
- // FIXME: do this locally
-
const Mesh& mesh(_C->mesh());
int q = mesh.topology().dim();
int n = _C->element().space_dimension();
@@ -285,7 +260,6 @@
solve(A, R_dT[e]->vector(), b, "cg");
}
end();
-
}
//-----------------------------------------------------------------------------
=== modified file 'dolfin/adaptivity/ErrorControl.h'
--- dolfin/adaptivity/ErrorControl.h 2010-11-04 16:16:16 +0000
+++ dolfin/adaptivity/ErrorControl.h 2010-11-05 12:06:31 +0000
@@ -1,7 +1,7 @@
// Copyright (C) 2010 Marie E. Rognes
//
// First added: 2010-08-19
-// Last changed: 2010-11-04
+// Last changed: 2010-11-05
#ifndef __ERROR_CONTROL_H
#define __ERROR_CONTROL_H
@@ -16,37 +16,112 @@
namespace dolfin
{
- // Forward declarations
- class Vector;
+ class BoundaryCondition;
class GoalFunctional;
- class BoundaryCondition;
class SpecialFacetFunction;
-
- class ErrorControl {
-
+ class Vector;
+
+ /// Generic Error Control class
+
+ class ErrorControl
+ {
public:
+ /// Create error control from variational problem and goal
+ /// functional. The sub-classed constructors must initialize all
+ /// forms and attach all relevant coefficients.
+ ///
+ /// *Arguments*
+ /// a
+ /// A (bilinear) _Form_.
+ /// L
+ /// A (linear) _Form_.
+ /// M
+ /// A _GoalFunctional_.
ErrorControl(const Form& a, const Form& L, const GoalFunctional& M);
+ /// Destructor.
~ErrorControl() { /* Do nothing */};
- // Estimate error
+ /// Estimate the error relative to the goal M of the discrete
+ /// approximation 'u' relative to the variational formulation by
+ /// evaluating the weak residual at an approximation to the dual
+ /// solution.
+ ///
+ /// *Arguments*
+ /// u
+ /// A _Function_.
+ /// bcs
+ /// A std::vector<const _BoundaryCondition_*>
double estimate_error(const Function& u,
std::vector<const BoundaryCondition*> bcs);
- // Compute error indicators
+ /// Compute error indicators
+ ///
+ /// *Arguments*
+ /// indicators
+ /// A _Vector_ (of length number of cells in mesh)
+ /// u
+ /// A _Function_
void compute_indicators(Vector& indicators, const Function& u);
- // Compute residual representation
+ /// Compute strong representation (strong cell and facet
+ /// residuals) of the weak residual.
+ ///
+ /// *Arguments*
+ /// R_T
+ /// A _Function_ (the strong cell residual)
+ /// R_dT
+ /// A _SpecialFacetFunction_ (the strong facet residual)
+ /// u
+ /// A _Function_
void residual_representation(Function& R_T,
SpecialFacetFunction& R_dT,
const Function& u);
+
+ /// Compute representation for the strong cell residual
+ /// from the weak residual
+ ///
+ /// *Arguments*
+ /// R_T
+ /// A _Function_ (the strong cell residual)
+ /// u
+ /// A _Function_
void compute_cell_residual(Function& R_T, const Function& u);
+
+ /// Compute representation for the strong facet residual from the
+ /// weak residual and the strong cell residual
+ ///
+ /// *Arguments*
+ /// R_dT
+ /// A _SpecialFacetFunction_ (the strong facet residual)
+ /// u
+ /// A _Function_
+ /// R_T
+ /// A _Function_ (the strong cell residual)
void compute_facet_residual(SpecialFacetFunction& R_dT,
const Function& u,
const Function& R_T);
+
+ /// Compute dual approximation defined by dual variational
+ /// problem and dual boundary conditions given by homogenized primal
+ /// boundary conditions.
+ ///
+ /// *Arguments*
+ /// z
+ /// A _Function_
+ /// bcs
+ /// A std::vector<const _BoundaryCondition_*>
void compute_dual(Function& z,
std::vector<const BoundaryCondition*> bcs);
+
+ /// Compute extrapolation with boundary conditions
+ ///
+ /// *Arguments*
+ /// z
+ /// A _Function_
+ /// bcs
+ /// A std::vector<const _BoundaryCondition_*>
void compute_extrapolation(const Function& z,
std::vector<const BoundaryCondition*> bcs);
@@ -62,10 +137,14 @@
// Functional for evaluating residual
boost::scoped_ptr<Form> _residual;
- // Bilinear and linear form for computing cell residual R_T
+ // Function (space)s for cell residual
+ // FIXME: meg is not sure if all of these (function spaces) are
+ // needed, just want to avoid things going out of scope.
boost::shared_ptr<FunctionSpace> _B;
boost::scoped_ptr<Function> _b_T;
- boost::shared_ptr<FunctionSpace> _DG_k; // This is
+ boost::shared_ptr<FunctionSpace> _DG_k;
+
+ // Bilinear and linear form for computing cell residual R_T
boost::scoped_ptr<Form> _a_R_T;
boost::scoped_ptr<Form> _L_R_T;
@@ -81,7 +160,7 @@
// Stored extrapolation
boost::shared_ptr<Function> _Ez_h;
- // References to input problem
+ // References to input variational problem
const Form& a;
const Form& L;
const GoalFunctional& M;
--- End Message ---
