I have an example of this situation based on step-12 which I have converted to perform explicit RK DG. It is based on step-39 also. If you look at this code
http://code.google.com/p/cfdlab/source/browse/trunk/deal.ii/step-12-RKDG/step-12.cc?spec=svn210&r=204 At line 234-238 I add my current solution vector to info_box. Then e.g., in function integrate_cell_term at line 326, I get "sol" out of cellinfo on line 330. This "sol" is the solution vector but already evaluated at the quadrature points of the current cell. Gradients can also be obtained similarly. praveen On Wed, Feb 9, 2011 at 9:58 AM, mihai alexe <[email protected]> wrote: > 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 > >
_______________________________________________ dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
