Mihai,
you are calling get_function_values explicitly in your cell assembling
routine, but this is not the intended use.
MeshWorker::IntegrationInfo can do this for you. Unfortunately, I have
not manage to get around to write an example program for this. You can
see how it works in the estimator of step 39.
Feel free to ask for more details if that does not help.
Best,
Guido
On 2/7/2011 2:57 PM, mihai alexe wrote:
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
_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
