CXX = /usr/local/gcc/gcc-2.8.1/bin/c++
endif
-cc-files = grid/grid_test.cc dof/dof_test.cc poisson/poisson.cc convergence/convergence.cc
+cc-files = grid/grid_test.cc dof/dof_test.cc poisson/poisson.cc \
+ convergence/convergence.cc error-estimation/error-estimation
o-files = $(cc-files:.cc=.o)
h-files = $(wildcard ../include/*.h)
@echo ================= Linking $@
@$(CXX) $(CXXFLAGS.g) -o $@ $< $(LIBS) ../../mia/control.o
+error-estimation/error-estimation: error-estimation/error-estimation.o $(LIBFILES.g)
+ @echo ================= Linking $@
+ @$(CXX) $(CXXFLAGS.g) -o $@ $< $(LIBS) ../../mia/control.o
-run: run_grid_test run_dof_test run_poisson_test run_convergence_test
+
+run: run_grid_test run_dof_test run_poisson_test run_convergence_test run_error_test
run_grid_test:
cd grid ; grid_test 4 ; mv grid.1 grid.4
cd convergence ; convergence
cd convergence ; gnuplot make_ps
+run_error_test:
+ cd error-estimation ; error-estimation
+ cd error-estimation ; gnuplot make_ps
+
clean:
cd grid ; rm -f grid.[1234] *.eps *.o *~ grid_test
cd dof ; rm -f grid.* sparsity.* *.o *~ dof_test
cd poisson ; make clean
cd convergence ; rm -f *.o convergence *~ *.eps gnuplot.*
+ cd error-estimation ; rm -f *.o error-estimation *~ *.eps gnuplot.*
-.PHONY: run run_grid_test run_dof_test run_poisson_test run_convergence_test clean
+.PHONY: run run_grid_test run_dof_test run_poisson_test run_convergence_test run_error_test clean
#Rule to generate the dependency file. This file is
n_dofs.push_back (dof->n_dofs());
cout << " Assembling matrices..." << endl;
- FEValues<dim>::UpdateStruct update_flags;
- update_flags.q_points = update_flags.gradients = true;
- update_flags.jacobians = update_flags.JxW_values = true;
+ UpdateFields update_flags = UpdateFields(update_q_points | update_gradients |
+ update_jacobians | update_JxW_values);
ProblemBase<dim>::DirichletBC dirichlet_bc;
dirichlet_bc[0] = boundary_values;
int main () {
PoissonProblem<2> problem;
- for (unsigned int level=1; level<10; ++level)
+ for (unsigned int level=1; level<9; ++level)
problem.run (level);
cout << endl << "Printing convergence history to <gnuplot.history>..." << endl;
cout << dof->n_dofs() << " degrees of freedom." << endl;
cout << " Assembling matrices..." << endl;
- FEValues<dim>::UpdateStruct update_flags;
- update_flags.q_points = update_flags.gradients = true;
- update_flags.jacobians = update_flags.JxW_values = true;
-
+ UpdateFields update_flags = UpdateFields(update_gradients | update_JxW_values);
ProblemBase<dim>::DirichletBC dirichlet_bc;
dirichlet_bc[0] = boundary_values;
assemble (equation, quadrature, fe, update_flags, dirichlet_bc);
template <int dim> class Quadrature;
+/**
+ Provide a set of flags which tells the #FEValues<>::reinit# function, which
+ fields are to be updated for each cell. E.g. if you do not need the
+ gradients since you want to assemble the mass matrix, you can switch that
+ off. By default, all flags are off, i.e. no reinitialization will be done.
+
+ A variable of this type has to be passed to the constructor of the
+ #FEValues# object. You can select more than one flag by concatenation
+ using the #|# (bitwise #or#) operator.
+ */
+enum UpdateFields {
+ /**
+ * Default: update nothing.
+ */
+ update_default = 0,
+ /**
+ * Compute quadrature points in real
+ * space (not on unit cell).
+ */
+ update_q_points = 1,
+ /**
+ * Transform gradients on unit cell to
+ * gradients on real cell.
+ */
+ update_gradients = 2,
+ /**
+ * Compute jacobian matrices of the
+ * transform between unit and real cell
+ * in the evaluation points.
+ */
+ update_jacobians = 4,
+ /**
+ * Compute the JxW values (Jacobian
+ * determinant at the quadrature point
+ * times the weight of this point).
+ */
+ update_JxW_values = 8,
+ /**
+ * Compute the points on the real cell
+ * on which the ansatz functions are
+ * located.
+ */
+ update_ansatz_points = 16
+};
+
+
+
+
/**
Represent a finite element evaluated with a specific quadrature rule.
This class is an optimization which avoids evaluating the shape functions
The unit cell is defined to be the tensor product of the interval $[0,1]$
in the present number of dimensions. In part of the literature, the convention
- is used that the unit cell be the tensor product of the intervale $[-1,1]$,
+ is used that the unit cell be the tensor product of the interval $[-1,1]$,
which is to distinguished properly.
Objects of this class store a multitude of different values needed to
The Jacobian matrix is defined to be
$$ J_{ij} = {d\xi_i \over dx_j} $$
- which is the form needed to compute the gradient on the real cell from
+ where the $\xi_i$ are the coordinates on the unit cell and the $x_i$ are
+ the coordinates on the real cell.
+ This is the form needed to compute the gradient on the real cell from
the gradient on the unit cell. If we want to transform the area element
$dx dy$ from the real to the unit cell, we have to take the determinant of
the inverse matrix, which is the reciprocal value of the determinant of the
The #FEValues# object keeps track of those fields which really need to
be computed, since the computation of the gradients of the ansatz functions
on each real cell can be quite an expensive thing if it is not needed. The
- object knows about which fields are needed by the #UpdateStruct# object
+ object knows about which fields are needed by the #UpdateFields# object
passed through the constructor. In debug mode, the accessor functions, which
return values from the different fields, check whether the required field
was initialized, thus avoiding use of unitialized data.
template <int dim>
class FEValues {
public:
- /**
- * Provide a structure which tells the
- * #reinit# function, which fields are
- * to be updated for each cell. E.g. if
- * you do not need the gradients since
- * you want to assemble the mass matrix,
- * you can switch that off. By default,
- * all flags are off, i.e. no
- * reinitialization will be done.
- *
- * A structure of this type has to be
- * passed to the constructor of the
- * #FEValues# object.
- */
- struct UpdateStruct {
- /**
- * Constructor. Sets all fields to
- * false.
- */
- UpdateStruct ();
- /**
- * Compute quadrature points in real
- * space (not on unit cell).
- */
- bool q_points;
- /**
- * Transform gradients on unit cell to
- * gradients on real cell.
- */
- bool gradients;
- /**
- * Compute jacobian matrices of the
- * transform between unit and real cell
- * in the evaluation points.
- */
- bool jacobians;
- /**
- * Compute the JxW values (Jacobian
- * determinant at the quadrature point
- * times the weight of this point).
- */
- bool JxW_values;
- /**
- * Compute the points on the real cell
- * on which the ansatz functions are
- * located.
- */
- bool ansatz_points;
- };
*/
FEValues (const FiniteElement<dim> &,
const Quadrature<dim> &,
- const UpdateStruct &);
+ const UpdateFields);
/**
* Return the value of the #i#th shape
* Store which fields are to be updated by
* the reinit function.
*/
- UpdateStruct update_flags;
+ UpdateFields update_flags;
};
inline
const vector<vector<Point<dim> > > &
FEValues<dim>::get_shape_grads () const {
- Assert (update_flags.gradients, ExcAccessToUninitializedField());
+ Assert (update_flags | update_gradients, ExcAccessToUninitializedField());
return shape_gradients;
};
inline
const vector<Point<dim> > &
FEValues<dim>::get_quadrature_points () const {
- Assert (update_flags.q_points, ExcAccessToUninitializedField());
+ Assert (update_flags | update_q_points, ExcAccessToUninitializedField());
return quadrature_points;
};
inline
const vector<Point<dim> > &
FEValues<dim>::get_ansatz_points () const {
- Assert (update_flags.ansatz_points, ExcAccessToUninitializedField());
+ Assert (update_flags | update_ansatz_points, ExcAccessToUninitializedField());
return ansatz_points;
};
inline
const vector<double> &
FEValues<dim>::get_JxW_values () const {
- Assert (update_flags.JxW_values, ExcAccessToUninitializedField());
+ Assert (update_flags | update_JxW_values, ExcAccessToUninitializedField());
return JxW_values;
};
dVector &rhs_vector,
const Quadrature<dim> &quadrature,
const FiniteElement<dim> &fe,
- const FEValues<dim>::UpdateStruct &update_flags);
+ const UpdateFields &update_flags);
/**
* Pointer to the dof handler object
* FEValues object need to be reinitialized
* on each cell.
*/
- const FEValues<dim>::UpdateStruct update_flags;
+ const UpdateFields update_flags;
};
* For what exactly happens here, refer to
* the general doc of this class.
*/
- virtual void assemble (const Equation<dim> &equation,
- const Quadrature<dim> &q,
- const FiniteElement<dim> &fe,
- const FEValues<dim>::UpdateStruct &update_flags,
- const DirichletBC &dirichlet_bc = DirichletBC());
+ virtual void assemble (const Equation<dim> &equation,
+ const Quadrature<dim> &q,
+ const FiniteElement<dim> &fe,
+ const UpdateFields &update_flags,
+ const DirichletBC &dirichlet_bc = DirichletBC());
/**
* Solve the system of equations.
-template <int dim>
-FEValues<dim>::UpdateStruct::UpdateStruct () :
- q_points(false),
- gradients(false),
- jacobians(false),
- JxW_values(false),
- ansatz_points(false) {};
-
-
-
/*------------------------------- FEValues -------------------------------*/
template <int dim>
FEValues<dim>::FEValues (const FiniteElement<dim> &fe,
const Quadrature<dim> &quadrature,
- const UpdateStruct &update_flags) :
+ const UpdateFields update_flags) :
n_quadrature_points(quadrature.n_quadrature_points),
total_dofs(fe.total_dofs),
shape_values(fe.total_dofs, quadrature.n_quadrature_points),
const unsigned int j) const {
Assert (i<(unsigned int)shape_values.m(), ExcInvalidIndex (i, shape_values.m()));
Assert (j<(unsigned int)shape_values.n(), ExcInvalidIndex (j, shape_values.n()));
- Assert (update_flags.gradients, ExcAccessToUninitializedField());
+ Assert (update_flags | update_gradients, ExcAccessToUninitializedField());
return shape_gradients[i][j];
};
template <int dim>
const Point<dim> & FEValues<dim>::quadrature_point (const unsigned int i) const {
Assert (i<n_quadrature_points, ExcInvalidIndex(i, n_quadrature_points));
- Assert (update_flags.q_points, ExcAccessToUninitializedField());
+ Assert (update_flags | update_q_points, ExcAccessToUninitializedField());
return quadrature_points[i];
};
template <int dim>
const Point<dim> & FEValues<dim>::ansatz_point (const unsigned int i) const {
Assert (i<ansatz_points.size(), ExcInvalidIndex(i, ansatz_points.size()));
- Assert (update_flags.ansatz_points, ExcAccessToUninitializedField());
+ Assert (update_flags | update_ansatz_points, ExcAccessToUninitializedField());
return ansatz_points[i];
};
template <int dim>
double FEValues<dim>::JxW (const unsigned int i) const {
Assert (i<n_quadrature_points, ExcInvalidIndex(i, n_quadrature_points));
- Assert (update_flags.JxW_values, ExcAccessToUninitializedField());
+ Assert (update_flags | update_JxW_values, ExcAccessToUninitializedField());
return JxW_values[i];
};
const FiniteElement<dim> &fe) {
// fill jacobi matrices and real
// quadrature points
- if (update_flags.jacobians || update_flags.q_points)
+ if ((update_flags | update_jacobians) ||
+ (update_flags | update_q_points))
fe.fill_fe_values (cell,
unit_quadrature_points,
jacobi_matrices,
- update_flags.jacobians,
+ update_flags | update_jacobians,
ansatz_points,
- update_flags.ansatz_points,
+ update_flags | update_ansatz_points,
quadrature_points,
- update_flags.q_points);
+ update_flags | update_q_points);
// compute gradients on real element if
// requested
- if (update_flags.gradients)
+ if (update_flags | update_gradients)
{
- Assert (update_flags.jacobians, ExcCannotInitializeField());
+ Assert (update_flags | update_jacobians, ExcCannotInitializeField());
for (unsigned int i=0; i<fe.total_dofs; ++i)
for (unsigned int j=0; j<n_quadrature_points; ++j)
// refer to the general doc for
// why we take the inverse of the
// determinant
- if (update_flags.JxW_values)
+ if (update_flags | update_JxW_values)
{
- Assert (update_flags.jacobians,
+ Assert (update_flags | update_jacobians,
ExcCannotInitializeField());
for (unsigned int i=0; i<n_quadrature_points; ++i)
JxW_values[i] = weights[i] / jacobi_matrices[i].determinant();
/*------------------------------- Explicit Instantiations -------------*/
-template struct FEValues<1>::UpdateStruct;
-template struct FEValues<2>::UpdateStruct;
-
template class FEValues<1>;
template class FEValues<2>;
dVector &rhs_vector,
const Quadrature<dim> &quadrature,
const FiniteElement<dim> &fe,
- const FEValues<dim>::UpdateStruct &update_flags) :
+ const UpdateFields &update_flags) :
dof(dof),
assemble_matrix(assemble_matrix),
assemble_rhs(assemble_rhs),
template <int dim>
-void ProblemBase<dim>::assemble (const Equation<dim> &equation,
- const Quadrature<dim> &quadrature,
- const FiniteElement<dim> &fe,
- const FEValues<dim>::UpdateStruct &update_flags,
- const DirichletBC &dirichlet_bc) {
+void ProblemBase<dim>::assemble (const Equation<dim> &equation,
+ const Quadrature<dim> &quadrature,
+ const FiniteElement<dim> &fe,
+ const UpdateFields &update_flags,
+ const DirichletBC &dirichlet_bc) {
Assert ((tria!=0) && (dof_handler!=0), ExcNoTriaSelected());
system_sparsity.reinit (dof_handler->n_dofs(),
difference.reinit (tria->n_active_cells());
- FEValues<dim>::UpdateStruct update_flags;
- update_flags.q_points = true;
- update_flags.jacobians = true;
- update_flags.JxW_values = true;
+ UpdateFields update_flags = UpdateFields (update_q_points |
+ update_jacobians |
+ update_JxW_values);
if ((norm==H1_seminorm) || (norm==H1_norm))
- update_flags.gradients = true;
+ update_flags = UpdateFields (update_flags | update_gradients);
FEValues<dim> fe_values(fe, q, update_flags);
// loop over all cells
CXX = /usr/local/gcc/gcc-2.8.1/bin/c++
endif
-cc-files = grid/grid_test.cc dof/dof_test.cc poisson/poisson.cc convergence/convergence.cc
+cc-files = grid/grid_test.cc dof/dof_test.cc poisson/poisson.cc \
+ convergence/convergence.cc error-estimation/error-estimation
o-files = $(cc-files:.cc=.o)
h-files = $(wildcard ../include/*.h)
@echo ================= Linking $@
@$(CXX) $(CXXFLAGS.g) -o $@ $< $(LIBS) ../../mia/control.o
+error-estimation/error-estimation: error-estimation/error-estimation.o $(LIBFILES.g)
+ @echo ================= Linking $@
+ @$(CXX) $(CXXFLAGS.g) -o $@ $< $(LIBS) ../../mia/control.o
-run: run_grid_test run_dof_test run_poisson_test run_convergence_test
+
+run: run_grid_test run_dof_test run_poisson_test run_convergence_test run_error_test
run_grid_test:
cd grid ; grid_test 4 ; mv grid.1 grid.4
cd convergence ; convergence
cd convergence ; gnuplot make_ps
+run_error_test:
+ cd error-estimation ; error-estimation
+ cd error-estimation ; gnuplot make_ps
+
clean:
cd grid ; rm -f grid.[1234] *.eps *.o *~ grid_test
cd dof ; rm -f grid.* sparsity.* *.o *~ dof_test
cd poisson ; make clean
cd convergence ; rm -f *.o convergence *~ *.eps gnuplot.*
+ cd error-estimation ; rm -f *.o error-estimation *~ *.eps gnuplot.*
-.PHONY: run run_grid_test run_dof_test run_poisson_test run_convergence_test clean
+.PHONY: run run_grid_test run_dof_test run_poisson_test run_convergence_test run_error_test clean
#Rule to generate the dependency file. This file is
n_dofs.push_back (dof->n_dofs());
cout << " Assembling matrices..." << endl;
- FEValues<dim>::UpdateStruct update_flags;
- update_flags.q_points = update_flags.gradients = true;
- update_flags.jacobians = update_flags.JxW_values = true;
+ UpdateFields update_flags = UpdateFields(update_q_points | update_gradients |
+ update_jacobians | update_JxW_values);
ProblemBase<dim>::DirichletBC dirichlet_bc;
dirichlet_bc[0] = boundary_values;
int main () {
PoissonProblem<2> problem;
- for (unsigned int level=1; level<10; ++level)
+ for (unsigned int level=1; level<9; ++level)
problem.run (level);
cout << endl << "Printing convergence history to <gnuplot.history>..." << endl;
cout << dof->n_dofs() << " degrees of freedom." << endl;
cout << " Assembling matrices..." << endl;
- FEValues<dim>::UpdateStruct update_flags;
- update_flags.q_points = update_flags.gradients = true;
- update_flags.jacobians = update_flags.JxW_values = true;
-
+ UpdateFields update_flags = UpdateFields(update_gradients | update_JxW_values);
ProblemBase<dim>::DirichletBC dirichlet_bc;
dirichlet_bc[0] = boundary_values;
assemble (equation, quadrature, fe, update_flags, dirichlet_bc);