Functional<number>::assemble(const DoFInfo<dim>& info)
{
for (unsigned int i=0;i<results.size();++i)
- results[i] += info.J[i];
+ results[i] += info.value(i);
}
{
for (unsigned int i=0;i<results.size();++i)
{
- results[i] += info1.J[i];
- results[i] += info2.J[i];
+ results[i] += info1.value(i);
+ results[i] += info2.value(i);
}
}
inline void
CellsAndFaces<number>::assemble(const DoFInfo<dim>& info)
{
- for (unsigned int i=0;i<info.J.size();++i)
+ for (unsigned int i=0;i<info.n_values();++i)
{
if (separate_faces &&
info.face_number != deal_II_numbers::invalid_unsigned_int)
- results(1)->block(i)(info.face->user_index()) += info.J[i];
+ results(1)->block(i)(info.face->user_index()) += info.value(i);
else
- results(0)->block(i)(info.cell->user_index()) += info.J[i];
+ results(0)->block(i)(info.cell->user_index()) += info.value(i);
}
}
CellsAndFaces<number>::assemble(const DoFInfo<dim>& info1,
const DoFInfo<dim>& info2)
{
- for (unsigned int i=0;i<info1.J.size();++i)
+ for (unsigned int i=0;i<info1.n_values();++i)
{
if (separate_faces)
{
- const double J = info1.J[i] + info2.J[i];
+ const double J = info1.value(i) + info2.value(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];
+ results(0)->block(i)(info1.cell->user_index()) += .5*info1.value(i);
+ results(0)->block(i)(info2.cell->user_index()) += .5*info2.value(i);
}
}
}
ResidualSimple<VECTOR>::assemble(const DoFInfo<dim>& info)
{
for (unsigned int k=0;k<residuals.size();++k)
- for (unsigned int i=0;i<info.R[k].block(0).size();++i)
- (*residuals(k))(info.indices[i]) += info.R[k].block(0)(i);
+ for (unsigned int i=0;i<info.vector(k).block(0).size();++i)
+ (*residuals(k))(info.indices[i]) += info.vector(k).block(0)(i);
}
{
for (unsigned int k=0;k<residuals.size();++k)
{
- for (unsigned int i=0;i<info1.R[k].block(0).size();++i)
- (*residuals(k))(info1.indices[i]) += info1.R[k].block(0)(i);
- for (unsigned int i=0;i<info2.R[k].block(0).size();++i)
- (*residuals(k))(info2.indices[i]) += info2.R[k].block(0)(i);
+ for (unsigned int i=0;i<info1.vector(k).block(0).size();++i)
+ (*residuals(k))(info1.indices[i]) += info1.vector(k).block(0)(i);
+ for (unsigned int i=0;i<info2.vector(k).block(0).size();++i)
+ (*residuals(k))(info2.indices[i]) += info2.vector(k).block(0)(i);
}
}
const DoFInfo<dim>& info)
{
for (unsigned int i=0;i<residuals.size();++i)
- assemble(*residuals(i), info.R[i], info.indices);
+ assemble(*residuals(i), info.vector(i), info.indices);
}
{
for (unsigned int i=0;i<residuals.size();++i)
{
- assemble(*residuals(i), info1.R[i], info1.indices);
- assemble(*residuals(i), info2.R[i], info2.indices);
+ assemble(*residuals(i), info1.vector(i), info1.indices);
+ assemble(*residuals(i), info2.vector(i), info2.indices);
}
}
inline void
MatrixSimple<MATRIX>::assemble(const DoFInfo<dim>& info)
{
- assemble(info.M1[0].matrix, info.indices, info.indices);
+ assemble(info.matrix(0,false).matrix, info.indices, info.indices);
}
MatrixSimple<MATRIX>::assemble(const DoFInfo<dim>& info1,
const DoFInfo<dim>& info2)
{
- assemble(info1.M1[0].matrix, info1.indices, info1.indices);
- assemble(info1.M2[0].matrix, info1.indices, info2.indices);
- assemble(info2.M1[0].matrix, info2.indices, info2.indices);
- assemble(info2.M2[0].matrix, info2.indices, info1.indices);
+ assemble(info1.matrix(0,false).matrix, info1.indices, info1.indices);
+ assemble(info1.matrix(0,true).matrix, info1.indices, info2.indices);
+ assemble(info2.matrix(0,false).matrix, info2.indices, info2.indices);
+ assemble(info2.matrix(0,true).matrix, info2.indices, info1.indices);
}
MGMatrixSimple<MATRIX>::assemble(const DoFInfo<dim>& info)
{
const unsigned int level = info.cell->level();
- assemble((*matrix)[level], info.M1[0].matrix, info.indices, info.indices);
+ assemble((*matrix)[level], info.matrix(0,false).matrix, info.indices, info.indices);
}
if (level1 == level2)
{
- assemble((*matrix)[level1], info1.M1[0].matrix, info1.indices, info1.indices);
- assemble((*matrix)[level1], info1.M2[0].matrix, info1.indices, info2.indices);
- assemble((*matrix)[level1], info2.M1[0].matrix, info2.indices, info2.indices);
- assemble((*matrix)[level1], info2.M2[0].matrix, info2.indices, info1.indices);
+ assemble((*matrix)[level1], info1.matrix(0,false).matrix, info1.indices, info1.indices);
+ assemble((*matrix)[level1], info1.matrix(0,true).matrix, info1.indices, info2.indices);
+ assemble((*matrix)[level1], info2.matrix(0,false).matrix, info2.indices, info2.indices);
+ assemble((*matrix)[level1], info2.matrix(0,true).matrix, info2.indices, info1.indices);
}
else
{
// Do not add info2.M1,
// which is done by
// the coarser cell
- assemble((*matrix)[level1], info1.M1[0].matrix, info1.indices, info1.indices);
- assemble_transpose((*flux_up)[level1],info1.M2[0].matrix, info2.indices, info1.indices);
- assemble((*flux_down)[level1], info2.M2[0].matrix, info2.indices, info1.indices);
+ assemble((*matrix)[level1], info1.matrix(0,false).matrix, info1.indices, info1.indices);
+ assemble_transpose((*flux_up)[level1],info1.matrix(0,true).matrix, info2.indices, info1.indices);
+ assemble((*flux_down)[level1], info2.matrix(0,true).matrix, info2.indices, info1.indices);
}
}
const unsigned int row = matrices[i]->row;
const unsigned int col = matrices[i]->column;
- assemble(matrices[i]->matrix, info.M1[i].matrix, row, col, info.indices, info.indices);
+ assemble(matrices[i]->matrix, info.matrix(i,false).matrix, row, col, info.indices, info.indices);
}
}
const unsigned int row = matrices[i]->row;
const unsigned int col = matrices[i]->column;
- assemble(matrices[i]->matrix, info1.M1[i].matrix, row, col, info1.indices, info1.indices);
- assemble(matrices[i]->matrix, info1.M2[i].matrix, row, col, info1.indices, info2.indices);
- assemble(matrices[i]->matrix, info2.M1[i].matrix, row, col, info2.indices, info2.indices);
- assemble(matrices[i]->matrix, info2.M2[i].matrix, row, col, info2.indices, info1.indices);
+ assemble(matrices[i]->matrix, info1.matrix(i,false).matrix, row, col, info1.indices, info1.indices);
+ assemble(matrices[i]->matrix, info1.matrix(i,true).matrix, row, col, info1.indices, info2.indices);
+ assemble(matrices[i]->matrix, info2.matrix(i,false).matrix, row, col, info2.indices, info2.indices);
+ assemble(matrices[i]->matrix, info2.matrix(i,true).matrix, row, col, info2.indices, info1.indices);
}
}
const unsigned int row = matrices[i]->row;
const unsigned int col = matrices[i]->column;
- assemble(matrices[i]->matrix[level], info.M1[i].matrix, row, col,
+ assemble(matrices[i]->matrix[level], info.matrix(i,false).matrix, row, col,
info.indices, info.indices, level, level);
}
}
if (level1 == level2)
{
- assemble(matrices[i]->matrix[level1], info1.M1[i].matrix, row, col, info1.indices, info1.indices, level1, level1);
- assemble(matrices[i]->matrix[level1], info1.M2[i].matrix, row, col, info1.indices, info2.indices, level1, level2);
- assemble(matrices[i]->matrix[level1], info2.M1[i].matrix, row, col, info2.indices, info2.indices, level2, level2);
- assemble(matrices[i]->matrix[level1], info2.M2[i].matrix, row, col, info2.indices, info1.indices, level2, level1);
+ assemble(matrices[i]->matrix[level1], info1.matrix(i,false).matrix, row, col, info1.indices, info1.indices, level1, level1);
+ assemble(matrices[i]->matrix[level1], info1.matrix(i,true).matrix, row, col, info1.indices, info2.indices, level1, level2);
+ assemble(matrices[i]->matrix[level1], info2.matrix(i,false).matrix, row, col, info2.indices, info2.indices, level2, level2);
+ assemble(matrices[i]->matrix[level1], info2.matrix(i,true).matrix, row, col, info2.indices, info1.indices, level2, level1);
}
else
{
// Do not add M22,
// which is done by
// the coarser cell
- assemble(matrices[i]->matrix[level1], info1.M1[i].matrix, row, col,
+ assemble(matrices[i]->matrix[level1], info1.matrix(i,false).matrix, row, col,
info1.indices, info1.indices, level1, level1);
- assemble(flux_up[i]->matrix[level1], info1.M2[i].matrix, row, col,
+ assemble(flux_up[i]->matrix[level1], info1.matrix(i,true).matrix, row, col,
info1.indices, info2.indices, level1, level2, true);
- assemble(flux_down[i]->matrix[level1], info2.M2[i].matrix, row, col,
+ assemble(flux_down[i]->matrix[level1], info2.matrix(i,true).matrix, row, col,
info2.indices, info1.indices, level2, level1);
}
}
*/
void reinit(const BlockIndices& local_sizes);
+ /**
+ * The number of scalar values.
+ */
+ unsigned int n_values () const;
+
+ /**
+ * The number of vectors.
+ */
+ unsigned int n_vectors () const;
+
+ /**
+ * The number of matrices.
+ */
+ unsigned int n_matrices () const;
+
+ /**
+ * Access scalar value at index
+ * @p i.
+ */
+ number& value(unsigned int i);
+
+ /**
+ * Read scalar value at index
+ * @p i.
+ */
+ number value(unsigned int i) const;
+
+ /**
+ * Access vector at index @p i.
+ */
+ BlockVector<number>& vector(unsigned int i);
+
+ /**
+ * Read vector at index @p i.
+ */
+ const BlockVector<number>& vector(unsigned int i) const;
+
+ /**
+ * Access matrix at index @p
+ * i. For results on internal
+ * faces, a true value for @p
+ * external refers to the flux
+ * between cells, while false
+ * refers to entries coupling
+ * inside the cell.
+ */
+ MatrixBlock<FullMatrix<number> >& matrix(unsigned int i, bool external = false);
+
+ /**
+ * Read matrix at index @p
+ * i. For results on internal
+ * faces, a true value for @p
+ * external refers to the flux
+ * between cells, while false
+ * refers to entries coupling
+ * inside the cell.
+ */
+ const MatrixBlock<FullMatrix<number> >& matrix(unsigned int i, bool external = false) const;
+
+ private:
/**
* The local numbers,
* computed on a cell or on a
/**
- * Basic info class only containing information on geometry and
- * degrees of freedom of the mesh object.
+ * A class containing information on geometry and degrees of freedom
+ * of a mesh object.
*
* The information in these objects is usually used by one of the
* Assembler classes. It is also the kind of information which is
*
* In addition to the information on degrees of freedom stored in this
* class, it also provides the local computation space for the worker
- * object operating on it. This space is provided by the base class
- * template DATATYPE. This base class will automatically
+ * object operating on it in LocalResults. This base class will automatically
* reinitialized on each cell, but initial setup is up to the user and
- * should be done when initialize() for this class is called. The
- * currently available base classes are
- * <ul>
- * <li> LocalVectors
- * <li> LocalMatrices
- * </ul>
+ * should be done when initialize() for this class is called.
*
* This class operates in two different modes, corresponding to the
* data models discussed in the Assembler namespace documentation.
* BlockInfo::initialize_local(). If this function has been used, or
* the vector has been changed from zero-length, then local dof
* indices stored in this object will automatically be renumbered to
- * reflect local block structure.
+ * reflect local block structure. This means, the first entries in
+ * #indices will refer to the first block of the system, then comes
+ * the second block and so on.
*
* The BlockInfo object is stored as a pointer. Therefore, if the
* block structure changes, for instance because of mesh refinement,
SmartPointer<const BlockInfo,DoFInfo<dim,spacedim> > block_info;
private:
- /// Fill index vector
+ /// Fill index vector with active indices
void get_indices(const typename DoFHandler<dim, spacedim>::cell_iterator& c);
/// Fill index vector with level indices
* @author Guido Kanschat, 2009
*/
template<int dim, class FEVALUESBASE, int spacedim = dim>
- class IntegrationInfo : public DoFInfo<dim, spacedim>
+ class IntegrationInfo
{
private:
/// vector of FEValues objects
std::vector<boost::shared_ptr<FEVALUESBASE> > fevalv;
public:
/**
- * Constructor forwarding
- * information to DoFInfo.
+ * Constructor.
*/
- IntegrationInfo(const BlockInfo& block_info);
-
- /**
- * Constructor forwarding
- * information to DoFInfo.
- */
- template <class DH>
- IntegrationInfo(const DH& dof_handler);
+ IntegrationInfo();
/**
* Build all internal
void initialize(const FiniteElement<dim,spacedim>& el,
const Mapping<dim,spacedim>& mapping,
const Quadrature<FEVALUES::integral_dimension>& quadrature,
- const UpdateFlags flags);
+ const UpdateFlags flags,
+ const BlockInfo* local_block_info = 0);
/**
* Initialize the data
* Reinitialize internal data
* structures for use on a cell.
*/
- template <class DHCellIterator>
- void reinit(const DHCellIterator& c);
-
+ void reinit(const DoFInfo<dim, spacedim>& i);
+
/**
- * Reinitialize internal data
- * structures for use on a face.
- */
- template <class DHCellIterator, class DHFaceIterator>
- void reinit(const DHCellIterator& c,
- const DHFaceIterator& f,
- const unsigned int fn);
-
- /**
- * Reinitialize internal data
- * structures for use on a subface.
- */
- template <class DHCellIterator, class DHFaceIterator>
- void reinit(const DHCellIterator& c,
- const DHFaceIterator& f,
- const unsigned int fn,
- const unsigned int sn);
-
-
- /**
- * @deprecated This is the
- * old version not using
- * VectorSelector.
- *
* Use the finite element
* functions in #global_data
* and fill the vectors
* #values, #gradients and
* #hessians.
*/
- void fill_local_data(const bool split_fevalues);
+ void fill_local_data(const DoFInfo<dim, spacedim>& info, const bool split_fevalues);
/**
* The global data vector
template <typename T>
IntegrationInfoBox(const T&);
- template <class WORKER, class ASSEMBLER>
+ template <class WORKER>
void initialize(const WORKER&,
- ASSEMBLER &assembler,
const FiniteElement<dim, spacedim>& el,
const Mapping<dim, spacedim>& mapping);
- template <class WORKER, class ASSEMBLER, typename VECTOR>
+ template <class WORKER, typename VECTOR>
void initialize(const WORKER&,
- ASSEMBLER &assembler,
const FiniteElement<dim, spacedim>& el,
const Mapping<dim, spacedim>& mapping,
const NamedData<VECTOR*>& data);
boost::shared_ptr<MeshWorker::VectorDataBase<dim, spacedim> > boundary_data;
boost::shared_ptr<MeshWorker::VectorDataBase<dim, spacedim> > face_data;
+ DoFInfo<dim, spacedim> dof_info;
CellInfo cell_info;
FaceInfo boundary_info;
FaceInfo face_info;
}
+ template <typename number>
+ inline
+ unsigned int
+ LocalResults<number>::n_values() const
+ {
+ return J.size();
+ }
+
+
+ template <typename number>
+ inline
+ unsigned int
+ LocalResults<number>::n_vectors() const
+ {
+ return R.size();
+ }
+
+
+ template <typename number>
+ inline
+ unsigned int
+ LocalResults<number>::n_matrices() const
+ {
+ return M1.size();
+ }
+
+
+ template <typename number>
+ inline
+ number&
+ LocalResults<number>::value(unsigned int i)
+ {
+ AssertIndexRange(i,J.size());
+ return J[i];
+ }
+
+
+ template <typename number>
+ inline
+ BlockVector<number>&
+ LocalResults<number>::vector(unsigned int i)
+ {
+ AssertIndexRange(i,R.size());
+ return R[i];
+ }
+
+
+ template <typename number>
+ inline
+ MatrixBlock<FullMatrix<number> >&
+ LocalResults<number>::matrix(unsigned int i, bool external)
+ {
+ if (external)
+ {
+ AssertIndexRange(i,M2.size());
+ return M2[i];
+ }
+ AssertIndexRange(i,M1.size());
+ return M1[i];
+ }
+
+ template <typename number>
+ inline
+ number
+ LocalResults<number>::value(unsigned int i) const
+ {
+ AssertIndexRange(i,J.size());
+ return J[i];
+ }
+
+
+ template <typename number>
+ inline
+ const BlockVector<number>&
+ LocalResults<number>::vector(unsigned int i) const
+ {
+ AssertIndexRange(i,R.size());
+ return R[i];
+ }
+
+
+ template <typename number>
+ inline
+ const MatrixBlock<FullMatrix<number> >&
+ LocalResults<number>::matrix(unsigned int i, bool external) const
+ {
+ if (external)
+ {
+ AssertIndexRange(i,M2.size());
+ return M2[i];
+ }
+ AssertIndexRange(i,M1.size());
+ return M1[i];
+ }
+
//----------------------------------------------------------------------//
template <int dim, int spacedim>
//----------------------------------------------------------------------//
- template <int dim, class FVB, int spacedim>
- template <class DH>
- IntegrationInfo<dim,FVB,spacedim>::IntegrationInfo(const DH& dof_handler)
- :
- DoFInfo<dim, spacedim>(dof_handler),
- fevalv(0),
- multigrid(false),
- global_data(boost::shared_ptr<VectorDataBase<dim, spacedim> >(new VectorDataBase<dim, spacedim>))
- {}
-
-
template <int dim, class FVB, int spacedim>
inline const FVB&
IntegrationInfo<dim,FVB,spacedim>::fe_values() const
template <int dim, class FVB, int spacedim>
- template <class DHCellIterator>
- inline void
- IntegrationInfo<dim,FVB,spacedim>::reinit(const DHCellIterator& c)
- {
- DoFInfo<dim,spacedim>::reinit(c);
- for (unsigned int i=0;i<fevalv.size();++i)
- {
- FVB& febase = *fevalv[i];
- FEValues<dim>& fe = dynamic_cast<FEValues<dim>&> (febase);
- fe.reinit(this->cell);
- }
-
- const bool split_fevalues = this->block_info != 0;
- fill_local_data(split_fevalues);
- }
-
-
- template <int dim, class FVB, int spacedim>
- template <class DHCellIterator, class DHFaceIterator>
inline void
- IntegrationInfo<dim,FVB,spacedim>::reinit(
- const DHCellIterator& c,
- const DHFaceIterator& f,
- const unsigned int fn)
+ IntegrationInfo<dim,FVB,spacedim>::reinit(const DoFInfo<dim, spacedim>& info)
{
- DoFInfo<dim,spacedim>::reinit(c, f, fn);
for (unsigned int i=0;i<fevalv.size();++i)
{
FVB& febase = *fevalv[i];
- FEFaceValues<dim>& fe = dynamic_cast<FEFaceValues<dim>&> (febase);
- fe.reinit(this->cell, fn);
+ if (info.sub_number != deal_II_numbers::invalid_unsigned_int)
+ {
+ // This is a subface
+ FESubfaceValues<dim>& fe = dynamic_cast<FESubfaceValues<dim>&> (febase);
+ fe.reinit(info.cell, info.face_number, info.sub_number);
+ }
+ else if (info.face_number != deal_II_numbers::invalid_unsigned_int)
+ {
+ // This is a face
+ FEFaceValues<dim>& fe = dynamic_cast<FEFaceValues<dim>&> (febase);
+ fe.reinit(info.cell, info.face_number);
+ }
+ else
+ {
+ // This is a cell
+ FEValues<dim>& fe = dynamic_cast<FEValues<dim>&> (febase);
+ fe.reinit(info.cell);
+ }
}
- const bool split_fevalues = this->block_info != 0;
- fill_local_data(split_fevalues);
+ const bool split_fevalues = info.block_info != 0;
+ fill_local_data(info, split_fevalues);
}
- template <int dim, class FVB, int spacedim>
- template <class DHCellIterator, class DHFaceIterator>
- inline void
- IntegrationInfo<dim,FVB,spacedim>::reinit(
- const DHCellIterator& c,
- const DHFaceIterator& f,
- const unsigned int fn,
- const unsigned int sn)
- {
- DoFInfo<dim,spacedim>::reinit(c, f, fn, sn);
- for (unsigned int i=0;i<fevalv.size();++i)
- {
- FVB& febase = *fevalv[i];
- FESubfaceValues<dim>& fe = dynamic_cast<FESubfaceValues<dim>&> (febase);
- fe.reinit(this->cell, fn, sn);
- }
-
- const bool split_fevalues = this->block_info != 0;
- fill_local_data(split_fevalues);
- }
-
//----------------------------------------------------------------------//
template <int dim, int sdim>
template <typename T>
IntegrationInfoBox<dim,sdim>::IntegrationInfoBox(const T& t)
:
- cell_info(t),
- boundary_info(t),
- face_info(t),
- subface_info(t),
- neighbor_info(t)
+ dof_info(t)
{}
template <int dim, int sdim>
- template <class WORKER, class ASSEMBLER>
+ template <class WORKER>
void
IntegrationInfoBox<dim,sdim>::
initialize(const WORKER& integrator,
- ASSEMBLER &assembler,
const FiniteElement<dim,sdim>& el,
const Mapping<dim,sdim>& mapping)
{
- assembler.initialize_info(cell_info, false);
- assembler.initialize_info(boundary_info, false);
- assembler.initialize_info(face_info, true);
- assembler.initialize_info(subface_info, true);
- assembler.initialize_info(neighbor_info, true);
-
cell_info.initialize<FEValues<dim,sdim> >(el, mapping, integrator.cell_quadrature,
integrator.cell_flags);
boundary_info.initialize<FEFaceValues<dim,sdim> >(el, mapping, integrator.boundary_quadrature,
template <int dim, int sdim>
- template <class WORKER, class ASSEMBLER, typename VECTOR>
+ template <class WORKER, typename VECTOR>
void
IntegrationInfoBox<dim,sdim>::initialize(
const WORKER& integrator,
- ASSEMBLER &assembler,
const FiniteElement<dim,sdim>& el,
const Mapping<dim,sdim>& mapping,
const NamedData<VECTOR*>& data)
{
- initialize(integrator, assembler, el, mapping);
+ initialize(integrator, el, mapping);
boost::shared_ptr<VectorData<VECTOR, dim, sdim> > p;
p = boost::shared_ptr<VectorData<VECTOR, dim, sdim> >(new VectorData<VECTOR, dim, sdim> (integrator.cell_selector));
//----------------------------------------------------------------------//
template<int dim, class FVB, int sdim>
- IntegrationInfo<dim,FVB,sdim>::IntegrationInfo(const BlockInfo& block_info)
+ IntegrationInfo<dim,FVB,sdim>::IntegrationInfo()
:
- DoFInfo<dim,sdim>(block_info),
fevalv(0),
multigrid(false),
global_data(boost::shared_ptr<VectorDataBase<dim, sdim> >(new VectorDataBase<dim, sdim>))
const FiniteElement<dim,sdim>& el,
const Mapping<dim,sdim>& mapping,
const Quadrature<FEVALUES::integral_dimension>& quadrature,
- const UpdateFlags flags)
+ const UpdateFlags flags,
+ const BlockInfo* block_info)
{
- if (this->block_info == 0 || this->block_info->local().size() == 0)
+ if (block_info == 0 || block_info->local().size() == 0)
{
fevalv.resize(1);
fevalv[0] = boost::shared_ptr<FVB> (
template<int dim, class FVB, int sdim>
void
- IntegrationInfo<dim,FVB,sdim>::fill_local_data(bool split_fevalues)
+ IntegrationInfo<dim,FVB,sdim>::fill_local_data(const DoFInfo<dim, sdim>& info, bool split_fevalues)
{
if (split_fevalues)
{
unsigned int comp = 0;
// Loop over all blocks
- for (unsigned int b=0;b<this->block_info->local().size();++b)
+ for (unsigned int b=0;b<info.block_info->local().size();++b)
{
- const unsigned int fe_no = this->block_info->base_element(b);
+ const unsigned int fe_no = info.block_info->base_element(b);
const FEValuesBase<dim>& fe = this->fe_values(fe_no);
const unsigned int n_comp = fe.get_fe().n_components();
- const unsigned int block_start = this->block_info->local().block_start(b);
- const unsigned int block_size = this->block_info->local().block_size(b);
+ const unsigned int block_start = info.block_info->local().block_start(b);
+ const unsigned int block_size = info.block_info->local().block_size(b);
- this->global_data->fill(values, gradients, hessians, fe, this->indices,
+ this->global_data->fill(values, gradients, hessians, fe, info.indices,
comp, n_comp, block_start, block_size);
comp += n_comp;
}
{
const FEValuesBase<dim>& fe = this->fe_values(0);
const unsigned int n_comp = fe.get_fe().n_components();
- this->global_data->fill(values, gradients, hessians, fe, this->indices,
- 0, n_comp, 0, this->indices.size());
+ this->global_data->fill(values, gradients, hessians, fe, info.indices,
+ 0, n_comp, 0, info.indices.size());
}
}
}
* @ingroup MeshWorker
* @author Guido Kanschat, 2009
*/
- template<class CELLINFO, class FACEINFO, class ITERATOR, typename ASSEMBLER>
+ template<class DOFINFO, class CELLINFO, class FACEINFO, class ITERATOR, typename ASSEMBLER>
void loop(ITERATOR begin,
typename identity<ITERATOR>::type end,
+ DOFINFO cell_dof_info,
CELLINFO& cell_info, FACEINFO& boundary_info,
FACEINFO& face_info, FACEINFO& subface_info, FACEINFO& neighbor_info,
- const std_cxx1x::function<void (CELLINFO &)> &cell_worker,
- const std_cxx1x::function<void (FACEINFO &)> &boundary_worker,
- const std_cxx1x::function<void (FACEINFO &, FACEINFO &)> &face_worker,
+ const std_cxx1x::function<void (DOFINFO&, CELLINFO &)> &cell_worker,
+ const std_cxx1x::function<void (DOFINFO&, FACEINFO &)> &boundary_worker,
+ const std_cxx1x::function<void (DOFINFO&, DOFINFO&, FACEINFO &, FACEINFO &)> &face_worker,
ASSEMBLER &assembler,
bool cells_first = true)
{
const bool integrate_boundary = (boundary_worker != 0);
const bool integrate_interior_face = (face_worker != 0);
+ ;
+ DOFINFO face_dof_info = cell_dof_info;
+ DOFINFO neighbor_dof_info = cell_dof_info;
+ assembler.initialize_info(cell_dof_info, false);
+ assembler.initialize_info(face_dof_info, true);
+ assembler.initialize_info(neighbor_dof_info, true);
+
// Loop over all cells
for (ITERATOR cell = begin; cell != end; ++cell)
{
// before faces
if (integrate_cell && cells_first)
{
- cell_info.reinit(cell);
- cell_worker(cell_info);
- assembler.assemble (cell_info);
+ cell_dof_info.reinit(cell);
+ cell_info.reinit(cell_dof_info);
+ cell_worker(cell_dof_info, cell_info);
+ assembler.assemble(cell_dof_info);
}
if (integrate_interior_face || integrate_boundary)
{
if (integrate_boundary)
{
- boundary_info.reinit(cell, face, face_no);
- boundary_worker(boundary_info);
- assembler.assemble (boundary_info);
+ cell_dof_info.reinit(cell, face, face_no);
+ boundary_info.reinit(cell_dof_info);
+ boundary_worker(cell_dof_info, boundary_info);
+ assembler.assemble(cell_dof_info);
}
}
else if (integrate_interior_face)
= cell->neighbor_of_coarser_neighbor(face_no);
typename ITERATOR::AccessorType::Container::face_iterator nface
= neighbor->face(neighbor_face_no.first);
- face_info.reinit(cell, face, face_no);
- subface_info.reinit(neighbor, nface, neighbor_face_no.first, neighbor_face_no.second);
+
+ face_dof_info.reinit(cell, face, face_no);
+ face_info.reinit(face_dof_info);
+ neighbor_dof_info.reinit(neighbor, nface,
+ neighbor_face_no.first, neighbor_face_no.second);
+ subface_info.reinit(neighbor_dof_info);
// Neighbor
// first to
// conform to
// old version
- face_worker(face_info, subface_info);
- assembler.assemble (face_info, subface_info);
+ face_worker(face_dof_info, neighbor_dof_info, face_info, subface_info);
+ assembler.assemble (face_dof_info, neighbor_dof_info);
}
else
{
unsigned int neighbor_face_no = cell->neighbor_of_neighbor(face_no);
Assert (neighbor->face(neighbor_face_no) == face, ExcInternalError());
// Regular interior face
- face_info.reinit(cell, face, face_no);
- neighbor_info.reinit(neighbor, neighbor->face(neighbor_face_no),
- neighbor_face_no);
-
- face_worker(face_info, neighbor_info);
- assembler.assemble (face_info, neighbor_info);
+ face_dof_info.reinit(cell, face, face_no);
+ face_info.reinit(face_dof_info);
+ neighbor_dof_info.reinit(neighbor, neighbor->face(neighbor_face_no),
+ neighbor_face_no);
+ neighbor_info.reinit(neighbor_dof_info);
+ face_worker(face_dof_info, neighbor_dof_info, face_info, neighbor_info);
+ assembler.assemble(face_dof_info, neighbor_dof_info);
}
}
} // faces
// have to be handled first
if (integrate_cell && !cells_first)
{
- cell_info.reinit(cell);
- cell_worker(cell_info);
- assembler.assemble (cell_info);
+ cell_dof_info.reinit(cell);
+ cell_info.reinit(cell_dof_info);
+ cell_worker(cell_dof_info, cell_info);
+ assembler.assemble (cell_dof_info);
}
}
}
void integration_loop(ITERATOR begin,
typename identity<ITERATOR>::type end,
IntegrationInfoBox<dim, dim>& box,
- const std_cxx1x::function<void (CELLINFO &)> &cell_worker,
- const std_cxx1x::function<void (FACEINFO &)> &boundary_worker,
- const std_cxx1x::function<void (FACEINFO &, FACEINFO &)> &face_worker,
+ const std_cxx1x::function<void (DoFInfo<dim>&, CELLINFO &)> &cell_worker,
+ const std_cxx1x::function<void (DoFInfo<dim>&, FACEINFO &)> &boundary_worker,
+ const std_cxx1x::function<void (DoFInfo<dim>&, DoFInfo<dim>&, FACEINFO &, FACEINFO &)> &face_worker,
ASSEMBLER &assembler,
bool cells_first = true)
{
- loop<CELLINFO,FACEINFO>
+ loop<DoFInfo<dim>,CELLINFO,FACEINFO>
(begin, end,
+ box.dof_info,
box.cell_info, box.boundary_info,
box.face_info,
box.subface_info, box.neighbor_info,
// $Id$
// Version: $Name$
//
-// Copyright (C) 2009 by the deal.II authors
+// Copyright (C) 2009, 2010 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
template void IntegrationInfo<deal_II_dimension, FEValuesBase<deal_II_dimension> >
::initialize<FEValues<deal_II_dimension> >(
const FiniteElement<deal_II_dimension>&, const Mapping<deal_II_dimension>&,
- const Quadrature<FEValues<deal_II_dimension>::integral_dimension>&, const UpdateFlags);
+ const Quadrature<FEValues<deal_II_dimension>::integral_dimension>&, const UpdateFlags, const BlockInfo*);
template void IntegrationInfo<deal_II_dimension, FEFaceValuesBase<deal_II_dimension> >
::initialize<FEFaceValues<deal_II_dimension> >(
const FiniteElement<deal_II_dimension>&, const Mapping<deal_II_dimension>&,
- const Quadrature<FEFaceValues<deal_II_dimension>::integral_dimension>&, const UpdateFlags);
+ const Quadrature<FEFaceValues<deal_II_dimension>::integral_dimension>&, const UpdateFlags, const BlockInfo*);
template void IntegrationInfo<deal_II_dimension, FEFaceValuesBase<deal_II_dimension> >
::initialize<FESubfaceValues<deal_II_dimension> >(
const FiniteElement<deal_II_dimension>&, const Mapping<deal_II_dimension>&,
- const Quadrature<FESubfaceValues<deal_II_dimension>::integral_dimension>&, const UpdateFlags);
+ const Quadrature<FESubfaceValues<deal_II_dimension>::integral_dimension>&, const UpdateFlags, const BlockInfo*);
}
#endif
// Here come the new include files
// for using the MeshWorker framework:
#include <numerics/mesh_worker.h>
+#include <numerics/mesh_worker_info.h>
#include <numerics/mesh_worker_loop.h>
#include <iostream>
// types of arguments, but have
// in fact other arguments
// already bound.
- static void integrate_cell_term (CellInfo& info);
- static void integrate_boundary_term (FaceInfo& info);
- static void integrate_face_term (FaceInfo& info1,
+ static void integrate_cell_term (MeshWorker::DoFInfo<dim>& dinfo, CellInfo& info);
+ static void integrate_boundary_term (MeshWorker::DoFInfo<dim>& dinfo, FaceInfo& info);
+ static void integrate_face_term (MeshWorker::DoFInfo<dim>& dinfo1,
+ MeshWorker::DoFInfo<dim>& dinfo2, FaceInfo& info1,
FaceInfo& info2);
};
// receives all the stuff we
// created so far.
MeshWorker::IntegrationInfoBox<dim> info_box(dof_handler);
- info_box.initialize(integration_worker, assembler, fe, mapping);
+ info_box.initialize(integration_worker, fe, mapping);
// Finally, the integration loop
// over all active cells
// (determined by the first
- // argument, which is an active iterator).
- MeshWorker::integration_loop<CellInfo,FaceInfo>
+ // argument, which is an active
+ // iterator).
+ MeshWorker::integration_loop<CellInfo, FaceInfo, dim>
(dof_handler.begin_active(), dof_handler.end(),
info_box,
&DGMethod<dim>::integrate_cell_term,
&DGMethod<dim>::integrate_boundary_term,
&DGMethod<dim>::integrate_face_term,
- assembler);
+ assembler, true);
}
// added soon).
template <int dim>
-void DGMethod<dim>::integrate_cell_term (CellInfo& info)
+void DGMethod<dim>::integrate_cell_term (MeshWorker::DoFInfo<dim>& dinfo, CellInfo& info)
{
// First, let us retrieve some of
// the objects used here from
// looks more complicated than
// might seem necessary.
const FEValuesBase<dim>& fe_v = info.fe_values();
- FullMatrix<double>& local_matrix = info.M1[0].matrix;
+ FullMatrix<double>& local_matrix = dinfo.matrix(0).matrix;
const std::vector<double> &JxW = fe_v.get_JxW_values ();
// With these objects, we continue
// FESubfaceValues, in order to get access to
// normal vectors.
template <int dim>
-void DGMethod<dim>::integrate_boundary_term (FaceInfo& info)
+void DGMethod<dim>::integrate_boundary_term (MeshWorker::DoFInfo<dim>& dinfo, FaceInfo& info)
{
const FEFaceValuesBase<dim>& fe_v = info.fe_values();
- FullMatrix<double>& local_matrix = info.M1[0].matrix;
- Vector<double>& local_vector = info.R[0].block(0);
+ FullMatrix<double>& local_matrix = dinfo.matrix(0).matrix;
+ Vector<double>& local_vector = dinfo.vector(0).block(0);
const std::vector<double> &JxW = fe_v.get_JxW_values ();
const std::vector<Point<dim> > &normals = fe_v.get_normal_vectors ();
// for each cell and two for coupling
// back and forth.
template <int dim>
-void DGMethod<dim>::integrate_face_term (FaceInfo& info1,
+void DGMethod<dim>::integrate_face_term (MeshWorker::DoFInfo<dim>& dinfo1,
+ MeshWorker::DoFInfo<dim>& dinfo2,
+ FaceInfo& info1,
FaceInfo& info2)
{
// For quadrature points, weights,
// interior couplings of that cell,
// while the second contains the
// couplings between cells.
- FullMatrix<double>& u1_v1_matrix = info1.M1[0].matrix;
- FullMatrix<double>& u2_v1_matrix = info1.M2[0].matrix;
- FullMatrix<double>& u1_v2_matrix = info2.M2[0].matrix;
- FullMatrix<double>& u2_v2_matrix = info2.M1[0].matrix;
+ FullMatrix<double>& u1_v1_matrix = dinfo1.matrix(0,false).matrix;
+ FullMatrix<double>& u2_v1_matrix = dinfo1.matrix(0,true).matrix;
+ FullMatrix<double>& u1_v2_matrix = dinfo2.matrix(0,true).matrix;
+ FullMatrix<double>& u2_v2_matrix = dinfo2.matrix(0,false).matrix;
// Here, following the previous
// functions, we would have the
class MatrixIntegrator : public Subscriptor
{
public:
- static void cell(typename MeshWorker::IntegrationWorker<dim>::CellInfo& info);
- static void bdry(typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info);
- static void face(typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info1,
+ static void cell(MeshWorker::DoFInfo<dim>& dinfo,
+ typename MeshWorker::IntegrationWorker<dim>::CellInfo& info);
+ static void bdry(MeshWorker::DoFInfo<dim>& dinfo,
+ typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info);
+ static void face(MeshWorker::DoFInfo<dim>& dinfo1,
+ MeshWorker::DoFInfo<dim>& dinfo2,
+ typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info1,
typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info2);
};
// degrees of freedom associated
// with the shape functions.
template <int dim>
-void MatrixIntegrator<dim>::cell(typename MeshWorker::IntegrationWorker<dim>::CellInfo& info)
+void MatrixIntegrator<dim>::cell(
+ MeshWorker::DoFInfo<dim>& dinfo,
+ typename MeshWorker::IntegrationWorker<dim>::CellInfo& info)
{
const FEValuesBase<dim>& fe = info.fe_values();
- FullMatrix<double>& local_matrix = info.M1[0].matrix;
+ FullMatrix<double>& local_matrix = dinfo.matrix(0,false).matrix;
for (unsigned int k=0; k<fe.n_quadrature_points; ++k)
for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
// On boundary faces, we use the
// Nitsche boundary condition
template <int dim>
-void MatrixIntegrator<dim>::bdry(typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info)
+void MatrixIntegrator<dim>::bdry(
+ MeshWorker::DoFInfo<dim>& dinfo,
+ typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info)
{
const FEFaceValuesBase<dim>& fe = info.fe_values();
- FullMatrix<double>& local_matrix = info.M1[0].matrix;
+ FullMatrix<double>& local_matrix = dinfo.matrix(0,false).matrix;
const unsigned int deg = fe.get_fe().tensor_degree();
- const double penalty = 2. * deg * (deg+1) * info.face->measure() / info.cell->measure();
+ const double penalty = 2. * deg * (deg+1) * dinfo.face->measure() / dinfo.cell->measure();
for (unsigned k=0;k<fe.n_quadrature_points;++k)
for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
template <int dim>
-void MatrixIntegrator<dim>::face(typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info1,
- typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info2)
+void MatrixIntegrator<dim>::face(
+ MeshWorker::DoFInfo<dim>& dinfo1,
+ MeshWorker::DoFInfo<dim>& dinfo2,
+ typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info1,
+ typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info2)
{
const FEFaceValuesBase<dim>& fe1 = info1.fe_values();
const FEFaceValuesBase<dim>& fe2 = info2.fe_values();
- FullMatrix<double>& matrix_v1u1 = info1.M1[0].matrix;
- FullMatrix<double>& matrix_v1u2 = info1.M2[0].matrix;
- FullMatrix<double>& matrix_v2u1 = info2.M2[0].matrix;
- FullMatrix<double>& matrix_v2u2 = info2.M1[0].matrix;
+ FullMatrix<double>& matrix_v1u1 = dinfo1.matrix(0,false).matrix;
+ FullMatrix<double>& matrix_v1u2 = dinfo1.matrix(0,true).matrix;
+ FullMatrix<double>& matrix_v2u1 = dinfo2.matrix(0,true).matrix;
+ FullMatrix<double>& matrix_v2u2 = dinfo2.matrix(0,false).matrix;
const unsigned int deg = fe1.get_fe().tensor_degree();
- double penalty1 = deg * (deg+1) * info1.face->measure() / info1.cell->measure();
- double penalty2 = deg * (deg+1) * info2.face->measure() / info2.cell->measure();
- if (info1.cell->has_children() ^ info2.cell->has_children())
+ double penalty1 = deg * (deg+1) * dinfo1.face->measure() / dinfo1.cell->measure();
+ double penalty2 = deg * (deg+1) * dinfo2.face->measure() / dinfo2.cell->measure();
+ if (dinfo1.cell->has_children() ^ dinfo2.cell->has_children())
{
- Assert (info1.face == info2.face, ExcInternalError());
- Assert (info1.face->has_children(), ExcInternalError());
-// Assert (info1.cell->has_children(), ExcInternalError());
+ Assert (dinfo1.face == dinfo2.face, ExcInternalError());
+ Assert (dinfo1.face->has_children(), ExcInternalError());
+// Assert (dinfo1.cell->has_children(), ExcInternalError());
penalty1 *= 2;
}
const double penalty = penalty1 + penalty2;
class RHSIntegrator : public Subscriptor
{
public:
- static void cell(typename MeshWorker::IntegrationWorker<dim>::CellInfo& info);
- static void bdry(typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info);
- static void face(typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info1,
+ static void cell(MeshWorker::DoFInfo<dim>& dinfo, typename MeshWorker::IntegrationWorker<dim>::CellInfo& info);
+ static void bdry(MeshWorker::DoFInfo<dim>& dinfo, typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info);
+ static void face(MeshWorker::DoFInfo<dim>& dinfo1,
+ MeshWorker::DoFInfo<dim>& dinfo2,
+ typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info1,
typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info2);
};
template <int dim>
-void RHSIntegrator<dim>::cell(typename MeshWorker::IntegrationWorker<dim>::CellInfo&)
+void RHSIntegrator<dim>::cell(MeshWorker::DoFInfo<dim>&, typename MeshWorker::IntegrationWorker<dim>::CellInfo&)
{}
template <int dim>
-void RHSIntegrator<dim>::bdry(typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info)
+void RHSIntegrator<dim>::bdry(MeshWorker::DoFInfo<dim>& dinfo, typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info)
{
const FEFaceValuesBase<dim>& fe = info.fe_values();
- Vector<double>& local_vector = info.R[0].block(0);
+ Vector<double>& local_vector = dinfo.vector(0).block(0);
std::vector<double> boundary_values(fe.n_quadrature_points);
exact_solution.value_list(fe.get_quadrature_points(), boundary_values);
const unsigned int deg = fe.get_fe().tensor_degree();
- const double penalty = 2. * deg * (deg+1) * info.face->measure() / info.cell->measure();
+ const double penalty = 2. * deg * (deg+1) * dinfo.face->measure() / dinfo.cell->measure();
for (unsigned k=0;k<fe.n_quadrature_points;++k)
for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
template <int dim>
-void RHSIntegrator<dim>::face(typename MeshWorker::IntegrationWorker<dim>::FaceInfo&,
+void RHSIntegrator<dim>::face(MeshWorker::DoFInfo<dim>&,
+ MeshWorker::DoFInfo<dim>&,
+ typename MeshWorker::IntegrationWorker<dim>::FaceInfo&,
typename MeshWorker::IntegrationWorker<dim>::FaceInfo&)
{}
class Estimator : public Subscriptor
{
public:
- static void cell(typename MeshWorker::IntegrationWorker<dim>::CellInfo& info);
- static void bdry(typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info);
- static void face(typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info1,
+ static void cell(MeshWorker::DoFInfo<dim>& dinfo, typename MeshWorker::IntegrationWorker<dim>::CellInfo& info);
+ static void bdry(MeshWorker::DoFInfo<dim>& dinfo, typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info);
+ static void face(MeshWorker::DoFInfo<dim>& dinfo1,
+ MeshWorker::DoFInfo<dim>& dinfo2,
+ typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info1,
typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info2);
};
template <int dim>
-void Estimator<dim>::cell(typename MeshWorker::IntegrationWorker<dim>::CellInfo& info)
+void Estimator<dim>::cell(MeshWorker::DoFInfo<dim>& dinfo, typename MeshWorker::IntegrationWorker<dim>::CellInfo& info)
{
const FEValuesBase<dim>& fe = info.fe_values();
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);
+ const double t = dinfo.cell->diameter() * trace(DDuh[k]);
+ dinfo.value(0) += t*t * fe.JxW(k);
}
}
template <int dim>
-void Estimator<dim>::bdry(typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info)
+void Estimator<dim>::bdry(MeshWorker::DoFInfo<dim>& dinfo, typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info)
{
const FEFaceValuesBase<dim>& fe = info.fe_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();
+ const double penalty = 2. * deg * (deg+1) * dinfo.face->measure() / dinfo.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])
+ dinfo.value(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,
+void Estimator<dim>::face(MeshWorker::DoFInfo<dim>& dinfo1,
+ MeshWorker::DoFInfo<dim>& dinfo2,
+ typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info1,
typename MeshWorker::IntegrationWorker<dim>::FaceInfo& info2)
{
const FEFaceValuesBase<dim>& fe = info1.fe_values();
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 penalty1 = deg * (deg+1) * dinfo1.face->measure() / dinfo1.cell->measure();
+ const double penalty2 = deg * (deg+1) * dinfo2.face->measure() / dinfo2.cell->measure();
const double penalty = penalty1 + penalty2;
- const double h = info1.face->measure();
+ const double h = dinfo1.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)
+ dinfo1.value(0) += (penalty * diff1*diff1 + h * diff2*diff2)
* fe.JxW(k);
}
- info2.J[0] = info1.J[0];
+ dinfo2.value(0) = dinfo1.value(0);
}
assembler.initialize(matrix);
MeshWorker::IntegrationInfoBox<dim> info_box(dof_handler);
- info_box.initialize(integration_worker, assembler, fe, mapping);
- MeshWorker::integration_loop<CellInfo, FaceInfo>(
+ info_box.initialize(integration_worker, fe, mapping);
+ MeshWorker::integration_loop<CellInfo, FaceInfo, dim>(
dof_handler.begin_active(), dof_handler.end(),
info_box,
&MatrixIntegrator<dim>::cell,
assembler.initialize(mg_matrix);
assembler.initialize_fluxes(mg_matrix_dg_up, mg_matrix_dg_down);
MeshWorker::IntegrationInfoBox<dim> info_box(mg_dof_handler);
- info_box.initialize(integration_worker, assembler, fe, mapping);
- MeshWorker::integration_loop<CellInfo, FaceInfo>(
+ info_box.initialize(integration_worker, fe, mapping);
+ MeshWorker::integration_loop<CellInfo, FaceInfo, dim>(
mg_dof_handler.begin(), mg_dof_handler.end(),
info_box,
&MatrixIntegrator<dim>::cell,
data.add(rhs, "RHS");
assembler.initialize(data);
MeshWorker::IntegrationInfoBox<dim> info_box(dof_handler);
- info_box.initialize(integration_worker, assembler, fe, mapping);
+ info_box.initialize(integration_worker, fe, mapping);
- MeshWorker::integration_loop<CellInfo, FaceInfo>(
+ MeshWorker::integration_loop<CellInfo, FaceInfo, dim>(
dof_handler.begin_active(), dof_handler.end(),
info_box,
&RHSIntegrator<dim>::cell,
out_data.add(est, "cells");
assembler.initialize(out_data, false);
MeshWorker::IntegrationInfoBox<dim> info_box(dof_handler);
- info_box.initialize(integration_worker, assembler, fe, mapping, solution_data);
- MeshWorker::integration_loop<CellInfo, FaceInfo> (
+ info_box.initialize(integration_worker, fe, mapping, solution_data);
+ MeshWorker::integration_loop<CellInfo, FaceInfo, dim> (
dof_handler.begin_active(), dof_handler.end(),
info_box,
&Estimator<dim>::cell,