Hello Guido,
Indeed, this is what I am doing. I did not know this was prohibited. Here is
what I am trying to do. I have a successive optimization problem. First, i
solve
the primal equation using the current guess for the control variables. Then, i
assemble and solve the dual, which depends on *both* the primal solution, and
the control variables. All variables are defined on the same mesh (and
associated with the same dof_handler object). How am I supposed to get the
primal and control variables (and gradients) at the quadrature points when
assembling my adjoint equation, if not through the
get_function_values/get_function_hessians routines? What I coded up looks as
follows:
template <int dim>
void NewtonOptimalitySystemSolver<dim>::integrate_cell_term_primal (DoFInfo&
dinfo, CellInfo& info)
{
const FEValuesBase<dim>& fe_v = info.fe_values();
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 ();
//use the FEValuesBase object to get the value of the optimal solution guess
//at the quadrature points...
std::vector<double> q_coeff_values(fe_v.n_quadrature_points);
//CRASHES HERE I GUESS...
fe_v.get_function_values(optimality_solution, q_coeff_values);
static PrimalRHSFunction<dim> f_rhs;
std::vector<double> f(fe_v.n_quadrature_points);
f_rhs.value_list(fe_v.get_quadrature_points(), q_coeff_values);
//assemble the cell term using the values
for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point)
for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j)
local_matrix(i,j) += q_coeff_values[point] *
fe_v.shape_grad(j,point) *
fe_v.shape_grad(i,point) *
JxW[point];
for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point)
for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
local_vector(i) += f[point] *
fe_v.shape_value(i,point) *
JxW[point];
}
I looked up step 39 and it seems we can only get the values and gradients for
one variable... but i need values/gradients/hessians for up to 3 different
discrete variables in my assembly routines. Is that possible without using
something like a FEFieldFunction object? (i am guessing that would introduce
large interpolation errors)...
thanks a lot for clearing things up!
-- Mihai
________________________________
Von: Guido Kanschat <[email protected]>
An: mihai alexe <[email protected]>
CC: deal.ii <[email protected]>
Gesendet: Dienstag, den 8. Februar 2011, 23:16:07 Uhr
Betreff: Re: [deal.II] problem with meshworker assembly
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
