* communicating the cell and
* face vectors.
*/
- typedef NamedData<SmartPointer<BlockVector<number>,CellsAndFaces<number> > > DataVectors;
+ typedef NamedData<BlockVector<number>*> DataVectors;
/**
* The initialization
*/
template <int dim>
void initialize_info(DoFInfo<dim>& info,
- bool interior_face);
+ bool interior_face) const;
/**
* Assemble the local values
};
+//----------------------------------------------------------------------//
+
+ template <typename number>
+ inline void
+ Functional<number>::initialize(unsigned int n)
+ {
+ results.resize(n);
+ }
+
+
+ template <typename number>
+ template<int dim>
+ inline void
+ Functional<number>::assemble(const DoFInfo<dim>& info)
+ {
+ for (unsigned int i=0;i<results.size();++i)
+ results[i] += info.J[i];
+ }
+
+
+ template <typename number>
+ template<int dim>
+ inline void
+ Functional<number>::assemble(const DoFInfo<dim>& info1,
+ const DoFInfo<dim>& info2)
+ {
+ for (unsigned int i=0;i<results.size();++i)
+ {
+ results[i] += info1.J[i];
+ results[i] += info2.J[i];
+ }
+ }
+
+
+ template <typename number>
+ inline number
+ Functional<number>::operator() (unsigned int i) const
+ {
+ AssertIndexRange(i, results.size());
+ return results[i];
+ }
+
+//----------------------------------------------------------------------//
+
+ template <typename number>
+ inline void
+ CellsAndFaces<number>::initialize(DataVectors& r, bool sep)
+ {
+ Assert(r.name(0) == "cells", typename DataVectors::ExcNameMismatch(0, "cells"));
+ if (sep)
+ {
+ Assert(r.name(1) == "faces", typename DataVectors::ExcNameMismatch(1, "faces"));
+ AssertDimension(r(0)->n_blocks(), r(1)->n_blocks());
+ }
+
+ results = r;
+ separate_faces = sep;
+ }
+
+
+ template <typename number>
+ template <int dim>
+ inline void
+ CellsAndFaces<number>::initialize_info(DoFInfo<dim>& info, bool) const
+ {
+ info.initialize_numbers(results(0)->n_blocks());
+ }
+
+
+ template <typename number>
+ template<int dim>
+ inline void
+ CellsAndFaces<number>::assemble(const DoFInfo<dim>& info)
+ {
+ for (unsigned int i=0;i<info.J.size();++i)
+ {
+ if (separate_faces &&
+ info.face_number != deal_II_numbers::invalid_unsigned_int)
+ results(1)->block(i)(info.face->user_index()) += info.J[i];
+ else
+ results(0)->block(i)(info.cell->user_index()) += info.J[i];
+ }
+ }
+
+
+ template <typename number>
+ template<int dim>
+ inline void
+ CellsAndFaces<number>::assemble(const DoFInfo<dim>& info1,
+ const DoFInfo<dim>& info2)
+ {
+ for (unsigned int i=0;i<info1.J.size();++i)
+ {
+ if (separate_faces)
+ {
+ const double J = info1.J[i] + info2.J[i];
+ results(1)->block(i)(info1.face->user_index()) += J;
+ if (info2.face != info1.face)
+ results(1)->block(i)(info2.face->user_index()) += J;
+ }
+ else
+ {
+ results(0)->block(i)(info1.cell->user_index()) += .5*info1.J[i];
+ results(0)->block(i)(info2.cell->user_index()) += .5*info2.J[i];
+ }
+ }
+ }
+
+
+
//----------------------------------------------------------------------//
template <class VECTOR>
MatrixBlock<FullMatrix<number> >& M,
unsigned int row, unsigned int col);
public:
- void initialize_vectors(const unsigned int n_vectors);
+ void initialize_numbers(const unsigned int n);
+ void initialize_vectors(const unsigned int n);
/**
* Allocate @p n local
* matrices. Additionally,
*/
template <typename T>
IntegrationInfoBox(const T&);
-
+
+ template <class VECTOR>
+ void initialize_data(const VECTOR*, const std::string& name,
+ bool values, bool gradients, bool hessians);
+
template <class WORKER>
void initialize(const WORKER&,
const FiniteElement<dim, spacedim>& el,
const Mapping<dim, spacedim>& mapping);
- template <class WORKER, class VECTOR,class P>
- void initialize(const WORKER&,
- const FiniteElement<dim, spacedim>& el,
- const Mapping<dim, spacedim>& mapping,
- const NamedData<SmartPointer<VECTOR,P> >& data);
// private:
-
+
+ boost::shared_ptr<MeshWorker::VectorDataBase<dim, spacedim> > data;
+
CellInfo cell_info;
FaceInfo bdry_info;
FaceInfo face_info;
//----------------------------------------------------------------------//
+ template <typename number>
+ inline void
+ LocalResults<number>::initialize_numbers(unsigned int n)
+ {
+ J.resize(n);
+ }
+
+
template <typename number>
inline void
LocalResults<number>::initialize_vectors(unsigned int n)
inline void
LocalResults<number>::reinit(const BlockIndices& bi)
{
+ for (unsigned int i=0;i<J.size();++i)
+ J[i] = 0.;
for (unsigned int i=0;i<R.size();++i)
R[i].reinit(bi);
for (unsigned int i=0;i<M1.size();++i)
subface_info(t),
neighbor_info(t)
{}
+
+
+ template <int dim, int sdim>
+ template <typename VECTOR>
+ void
+ IntegrationInfoBox<dim,sdim>::initialize_data(
+ const VECTOR* v, const std::string& name,
+ bool values, bool gradients, bool hessians)
+ {
+ boost::shared_ptr<VectorData<VECTOR, dim, sdim> >
+ p = boost::shared_ptr<VectorData<VECTOR, dim, sdim> >(new VectorData<VECTOR, dim, sdim> ());
+ p->add(name, values, gradients, hessians);
+ p->initialize(v, name);
+ data = p;
+ cell_info.initialize_data(data);
+ bdry_info.initialize_data(data);
+ face_info.initialize_data(data);
+ subface_info.initialize_data(data);
+ neighbor_info.initialize_data(data);
+ }
+
+
template <int dim, int sdim>
template <class WORKER>
void
cell_info.initialize<FEValues<dim,sdim> >(el, mapping, integrator.cell_quadrature,
integrator.cell_flags);
bdry_info.initialize<FEFaceValues<dim,sdim> >(el, mapping, integrator.bdry_quadrature,
- integrator.face_flags);
+ integrator.bdry_flags);
face_info.initialize<FEFaceValues<dim,sdim> >(el, mapping, integrator.face_quadrature,
integrator.face_flags);
subface_info.initialize<FESubfaceValues<dim,sdim> >(el, mapping, integrator.face_quadrature,
integrator.face_flags);
neighbor_info.initialize<FEFaceValues<dim,sdim> >(el, mapping, integrator.face_quadrature,
integrator.ngbr_flags);
- }
-
-
- template <int dim, int sdim>
- template <class WORKER, class VECTOR, class P>
- void
- IntegrationInfoBox<dim,sdim>::initialize(
- const WORKER& integrator,
- const FiniteElement<dim,sdim>& el,
- const Mapping<dim,sdim>& mapping,
- const NamedData<SmartPointer<VECTOR,P> >& data)
- {
- cell_info.initialize_data(data);
- bdry_info.initialize_data(data);
- face_info.initialize_data(data);
- subface_info.initialize_data(data);
- neighbor_info.initialize_data(data);
-
- initialize(integrator, el, mapping);
- }
+ }
}
// For all components
for (unsigned int j=0;j<values[i].size();++j)
{
- values[i][j].resize(el.n_components());
+ values[i][j].resize(quadrature.size());
}
}
// For all components
for (unsigned int j=0;j<gradients[i].size();++j)
{
- gradients[i][j].resize(el.n_components());
+ gradients[i][j].resize(quadrature.size());
}
}
// For all components
for (unsigned int j=0;j<hessians[i].size();++j)
{
- hessians[i][j].resize(el.n_components());
+ hessians[i][j].resize(quadrature.size());
}
}
}
template <class STREAM, typename DATA>
void print (STREAM& s, const NamedData<DATA>& v) const;
+ /**
+ * Print the number of selections to the stream.
+ */
+ template <class STREAM>
+ void print (STREAM& s) const;
+
protected:
/**
* Selection of the vectors
unsigned int size) const;
};
+
/**
* Based on VectorSelector, this is the class that implements the
- * function VectorDataBase::fill() for a certain type of vector.
+ * function VectorDataBase::fill() for a certain type of vector, using
+ * NamedData to identify vectors by name.
*
* @author Guido Kanschat, 2009
*/
{
public:
void initialize(const NamedData<VECTOR*>&);
+ void initialize(const VECTOR*, const std::string& name);
virtual void fill(
std::vector<std::vector<std::vector<double> > >& values,
unsigned int start,
unsigned int size) const;
private:
- NamedData<SmartPointer<VECTOR,VectorData<VECTOR,dim,spacedim> > > data;
+ NamedData<SmartPointer<const VECTOR,VectorData<VECTOR,dim,spacedim> > > data;
};
+
//----------------------------------------------------------------------//
}
+ template <class STREAM>
+ inline void
+ VectorSelector::print(STREAM& s) const
+ {
+ s << "values: " << n_values()
+ << " gradients: " << n_gradients()
+ << " hessians: " << n_hessians()
+ << std::endl;
+ }
+
+
template <class STREAM, typename DATA>
inline void
VectorSelector::print(STREAM& s, const NamedData<DATA>& v) const
}
+ template <class VECTOR, int dim, int spacedim>
+ void
+ VectorData<VECTOR, dim, spacedim>::initialize(const VECTOR* v, const std::string& name)
+ {
+ SmartPointer<const VECTOR,VectorData<VECTOR, dim, spacedim> > p = v;
+ data.add(p, name);
+ VectorSelector::initialize(data);
+ }
+
+
template <class VECTOR, int dim, int spacedim>
void
VectorData<VECTOR, dim, spacedim>::fill(
namespace Assembler
{
//----------------------------------------------------------------------//
-
- template <typename number>
- inline void
- Functional<number>::initialize(unsigned int n)
- {
- results.resize(n);
- }
-
-
- template <typename number>
- template<int dim>
- inline void
- Functional<number>::assemble(const DoFInfo<dim>& info)
- {
- for (unsigned int i=0;i<results.size();++i)
- results[i] += info.J[i];
- }
-
- template <typename number>
- template<int dim>
- inline void
- Functional<number>::assemble(const DoFInfo<dim>& info1,
- const DoFInfo<dim>& info2)
- {
- for (unsigned int i=0;i<results.size();++i)
- {
- results[i] += info1.J[i];
- results[i] += info2.J[i];
- }
- }
-
-
- template <typename number>
- inline number
- Functional<number>::operator() (unsigned int i) const
- {
- AssertIndexRange(i, results.size());
- return results[i];
- }
-
//----------------------------------------------------------------------//
- template <typename number>
- inline void
- CellsAndFaces<number>::initialize(DataVectors& r, bool sep)
- {
- Assert(r.name(0) == "cells", DataVectors::ExcNameMismatch(0, "cells"));
- if (sep)
- Assert(r.name(1) == "faces", DataVectors::ExcNameMismatch(1, "faces"));
- AssertDimension(r(0).n_blocks(), r(1).n_blocks());
-
- results = r;
- separate_faces = sep;
- }
-
-
- template <typename number>
- template<int dim>
- inline void
- CellsAndFaces<number>::assemble(const DoFInfo<dim>& info)
- {
- for (unsigned int i=0;i<info.J.size();++i)
- {
- if (separate_faces &&
- info.face_number != deal_II_numbers::invalid_unsigned_int)
- results.vector(1).block(i)(info.face->user_index()) += info.J[i];
- else
- results.vector(0).block(i)(info.cell->user_index()) += info.J[i];
- }
- }
-
-
- template <typename number>
- template<int dim>
- inline void
- CellsAndFaces<number>::assemble(const DoFInfo<dim>& info1,
- const DoFInfo<dim>& info2)
- {
- for (unsigned int i=0;i<info1.J.size();++i)
- {
- if (separate_faces)
- {
- const double J = info1.J[i] + info2.J[i];
- results.vector(1).block(i)(info1.face->user_index()) += J;
- if (info2.face != info1.face)
- results.vector(1).block(i)(info2.face->user_index()) += J;
- }
- else
- {
- results.vector(0).block(i)(info1.cell->user_index()) += info1.J[i];
- results.vector(0).block(i)(info2.cell->user_index()) += info2.J[i];
- }
- }
- }
-
-
-//----------------------------------------------------------------------//
-
-
-
}
}
#include <lac/solver_cg.h>
#include <lac/precondition.h>
#include <lac/precondition_block.h>
+#include <lac/block_vector.h>
// Include files for setting up the
// mesh
#include <grid/grid_generator.h>
+#include <grid/grid_refinement.h>
// Include files for FiniteElement
// classes and DoFHandler.
};
+ // On each cell, we integrate the
+ // Dirichlet form. All local
+ // integrations consist of nested
+ // loops, first over all quadrature
+ // points and the iner loops over the
+ // degrees of freedom associated
+ // with the shape functions.
template <int dim>
void MatrixIntegrator<dim>::cell(typename MeshWorker::IntegrationWorker<dim>::CellInfo& info) const
{
* fe.JxW(k);
}
-
+ // On boundary faces, we use the
+ // Nitsche boundary condition
template <int dim>
void MatrixIntegrator<dim>::bdry(typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info) const
{
{
Assert (info1.face == info2.face, ExcInternalError());
Assert (info1.face->has_children(), ExcInternalError());
- Assert (info1.cell->has_children(), ExcInternalError());
+// Assert (info1.cell->has_children(), ExcInternalError());
penalty1 *= 2;
}
const double penalty = penalty1 + penalty2;
{}
+template <int dim>
+class Estimator : public Subscriptor
+{
+ public:
+ void cell(typename MeshWorker::IntegrationWorker<dim>::CellInfo& info) const;
+ void bdry(typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info) const;
+ void face(typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info1,
+ typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info2) const;
+};
+
+
+template <int dim>
+void Estimator<dim>::cell(typename MeshWorker::IntegrationWorker<dim>::CellInfo& info) const
+{
+ const FEValuesBase<dim>& fe = info.fe();
+
+ const std::vector<Tensor<2,dim> >& DDuh = info.hessians[0][0];
+ for (unsigned k=0;k<fe.n_quadrature_points;++k)
+ {
+ const double t = info.cell->diameter() * trace(DDuh[k]);
+ info.J[0] += t*t * fe.JxW(k);
+ }
+}
+
+
+template <int dim>
+void Estimator<dim>::bdry(typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info) const
+{
+ const FEFaceValuesBase<dim>& fe = info.fe();
+
+ std::vector<double> boundary_values(fe.n_quadrature_points);
+ exact_solution.value_list(fe.get_quadrature_points(), boundary_values);
+
+ const std::vector<double>& uh = info.values[0][0];
+
+ const unsigned int deg = fe.get_fe().tensor_degree();
+ const double penalty = 2. * deg * (deg+1) * info.face->measure() / info.cell->measure();
+
+ for (unsigned k=0;k<fe.n_quadrature_points;++k)
+ info.J[0] += penalty * (boundary_values[k] - uh[k]) * (boundary_values[k] - uh[k])
+ * fe.JxW(k);
+}
+
+
+template <int dim>
+void Estimator<dim>::face(typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info1,
+ typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info2) const
+{
+ const FEFaceValuesBase<dim>& fe = info1.fe();
+ const std::vector<double>& uh1 = info1.values[0][0];
+ const std::vector<double>& uh2 = info2.values[0][0];
+ const std::vector<Tensor<1,dim> >& Duh1 = info1.gradients[0][0];
+ const std::vector<Tensor<1,dim> >& Duh2 = info2.gradients[0][0];
+
+ const unsigned int deg = fe.get_fe().tensor_degree();
+ const double penalty1 = deg * (deg+1) * info1.face->measure() / info1.cell->measure();
+ const double penalty2 = deg * (deg+1) * info2.face->measure() / info2.cell->measure();
+ const double penalty = penalty1 + penalty2;
+ const double h = info1.face->measure();
+
+ for (unsigned k=0;k<fe.n_quadrature_points;++k)
+ {
+ double diff1 = uh1[k] - uh2[k];
+ double diff2 = fe.normal_vector(k) * Duh1[k] - fe.normal_vector(k) * Duh2[k];
+ info1.J[0] += (penalty * diff1*diff1 + h * diff2*diff2)
+ * fe.JxW(k);
+ }
+ info2.J[0] = info1.J[0];
+}
+
+
// @sect3{The main class}
template <int dim>
class Step39
void assemble_mg_matrix ();
void assemble_right_hand_side ();
void error ();
- void estimate ();
+ double estimate ();
void solve ();
void refine_grid ();
void output_results (const unsigned int cycle) const;
SparseMatrix<double> matrix;
Vector<double> solution;
Vector<double> right_hand_side;
-
+ BlockVector<double> estimates;
+
MGLevelObject<SparsityPattern> mg_sparsity;
MGLevelObject<SparsityPattern> mg_sparsity_dg_interface;
MGLevelObject<SparseMatrix<double> > mg_matrix;
:
fe(fe),
mg_dof_handler(triangulation),
- dof_handler(mg_dof_handler)
+ dof_handler(mg_dof_handler),
+ estimates(1)
{
GridGenerator::hyper_L(triangulation, -1, 1);
- triangulation.begin()->set_refine_flag();
- triangulation.execute_coarsening_and_refinement();
}
QGauss<dim> quadrature(n_gauss_points);
VectorTools::integrate_difference(mapping, dof_handler, solution, exact_solution,
cell_errors, quadrature, VectorTools::L2_norm);
- deallog << "Error " << cell_errors.l2_norm() << std::endl;
+ deallog << "Error " << cell_errors.l2_norm() << std::endl;
+}
+
+
+template <int dim>
+double
+Step39<dim>::estimate()
+{
+ estimates.block(0).reinit(triangulation.n_active_cells());
+ unsigned int i=0;
+ for (typename Triangulation<dim>::active_cell_iterator cell = triangulation.begin_active();
+ cell != triangulation.end();++cell,++i)
+ cell->set_user_index(i);
+
+ const Estimator<dim> local;
+ MeshWorker::AssemblingIntegrator<dim, MeshWorker::Assembler::CellsAndFaces<double>, Estimator<dim> >
+ integrator(local);
+ const unsigned int n_gauss_points = dof_handler.get_fe().tensor_degree()+1;
+ integrator.initialize_gauss_quadrature(n_gauss_points, n_gauss_points+1, n_gauss_points);
+ UpdateFlags update_flags = update_values | update_gradients;
+ integrator.add_update_flags(update_flags | update_hessians, true, true, true, true);
+// integrator.add_update_flags(update_hessians, true, false, false, false);
+ integrator.add_update_flags(update_quadrature_points, false, true, false, false);
+
+ NamedData<BlockVector<double>* > out_data;
+ BlockVector<double>* est = &estimates;
+ out_data.add(est, "cells");
+ integrator.initialize(out_data, false);
+ MeshWorker::IntegrationInfoBox<dim> info_box(dof_handler);
+ info_box.initialize_data(&solution, std::string("solution"), true, true, true);
+ info_box.initialize(integrator, fe, mapping);
+ MeshWorker::integration_loop(dof_handler.begin_active(), dof_handler.end(), info_box, integrator);
+ return estimates.block(0).l2_norm();
}
{
// Output of the solution in
// gnuplot format.
- std::string filename = "sol-";
- filename += ('0' + cycle);
- Assert (cycle < 10, ExcInternalError());
-
+ char * fn = new char[100];
+ sprintf(fn, "sol-%02d", cycle);
+
+ std::string filename(fn);
filename += ".gnuplot";
std::cout << "Writing solution to <" << filename << ">..."
<< std::endl << std::endl;
DataOut<dim> data_out;
data_out.attach_dof_handler (dof_handler);
data_out.add_data_vector (solution, "u");
+ data_out.add_data_vector (estimates.block(0), "est");
data_out.build_patches ();
{
deallog << "Step " << s << std::endl;
if (s != 0)
- triangulation.refine_global(1);
+ {
+ if (estimates.block(0).size() == 0)
+ triangulation.refine_global(1);
+ else
+ {
+ GridRefinement::refine_and_coarsen_fixed_fraction (triangulation,
+ estimates.block(0),
+ 0.5, 0.0);
+ triangulation.execute_coarsening_and_refinement ();
+ }
+ }
deallog << "Triangulation "
<< triangulation.n_active_cells() << " cells, "
assemble_right_hand_side();
deallog << "Solve" << std::endl;
solve();
- deallog << "Error" << std::endl;
error();
+ deallog << "Estimate " << estimate() << std::endl;
output_results(s);
}
}
int main()
{
- FE_DGQ<2> fe1(1);
+ FE_DGQ<2> fe1(3);
Step39<2> test1(fe1);
- test1.run(7);
+ test1.run(13);
}