Hello,
I am running into some problems when assembling a matrix with the MeshWorker
framework. My code looks as follows:
//setup part
CompressedSparsityPattern c_sparsity(dof_handler.n_dofs());
DoFTools::make_flux_sparsity_pattern (dof_handler, c_sparsity);
sparsity_pattern.copy_from(c_sparsity);
uuv_system_matrix.reinit (sparsity_pattern);
primal_solution.reinit (dof_handler.n_dofs());
primal_rhs.reinit (dof_handler.n_dofs());
//assembly part:
MeshWorker::IntegrationInfoBox<dim> info_box;
const unsigned int n_gauss_points = dof_handler.get_fe().degree+1;
info_box.initialize_gauss_quadrature(n_gauss_points,
n_gauss_points,
n_gauss_points);
info_box.initialize_update_flags();
UpdateFlags update_flags = update_quadrature_points |
update_values |
update_gradients;
info_box.add_update_flags(update_flags, true, true, true, true);
info_box.initialize(fe, mapping);
MeshWorker::DoFInfo<dim> dof_info(dof_handler);
MeshWorker::Assembler::SystemSimple<SparseMatrix<double>, Vector<double> >
assembler;
assembler.initialize(system_matrix, primal_rhs);
MeshWorker::integration_loop<dim, dim>
(dof_handler.begin_active(), dof_handler.end(),
dof_info, info_box,
boost::bind(&NewtonOptimalitySystemSolver<dim>::integrate_cell_term_primal,
this, _1, _2),
boost::bind(&NewtonOptimalitySystemSolver<dim>::integrate_boundary_term_primal,
this, _1, _2),
boost::bind(&NewtonOptimalitySystemSolver<dim>::integrate_face_term_primal,
this, _1, _2, _3, _4),
assembler, true);
where the 3 functions assemble each part of the weak DG formulation (cell,
boundary, and interior faces). When I run the code, I get the following error:
--------------------------------------------------------
An error occurred in line <1585> of file
</home/mihai/Software/deal.II.7/deal.II/source/fe/fe_values.cc> in function
unsigned int dealii::FEValuesBase<dim,
spacedim>::TriaCellIterator::n_dofs_for_dof_handler() const [with int dim = 2,
int spacedim = 2]
The violated condition was:
false
The name and call sequence of the exception was:
ExcMessage (message_string)
Additional Information:
You have previously called the FEValues::reinit function with a
cell iterator of type Triangulation<dim,spacedim>::cell_iterator. However,
when you do this, you cannot call some functions in the FEValues
class, such as the get_function_values/gradients/hessians
functions. If you need these functions, then you need to call
FEValues::reinit with an iterator type that allows to extract
degrees of freedom, such as DoFHandler<dim,spacedim>::cell_iterator.
Stacktrace:
-----------
#0 /home/mihai/Software/deal.II.7/deal.II/lib/libdeal_II.g.so.7.0.0:
dealii::FEValuesBase<2, 2>::TriaCellIterator::n_dofs_for_dof_handler() const
#1 /home/mihai/Software/deal.II.7/deal.II/lib/libdeal_II.g.so.7.0.0: void
dealii::FEValuesBase<2, 2>::get_function_values<dealii::Vector<double>,
double>(dealii::Vector<double> const&, std::vector<double,
std::allocator<double> >&) const
#2 ./LbfgsbOptimizer:
NewtonOptimalitySystemSolver<2>::integrate_cell_term_primal(dealii::MeshWorker::DoFInfo<2,
2, double>&, dealii::MeshWorker::IntegrationInfo<2, 2>&)
... (CONTINUES)
The error message is pretty clear to me, but I am not calling FEValues::reinit
explicitly anywhere in my code - that's what is confusing. Where could this
call
come from? Before setting up the system, I just create the triangulation,
distribute the DoFs, and project a known function onto my mesh (with
VectorTools::project) - that's it. What could be causing this problem?
-- thanks!
-- Mihai
_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
