I'm not sure if this is helpful but as a diagnostic, replacing
MeshWorker::loop with this loop that I wrote myself:
const QGauss<dim> quadrature_formula(fe.degree + 1);
const unsigned int n_q_points = quadrature_formula.size();
FEValues<dim> fe_values(fe, quadrature_formula, update_values |
update_JxW_values);
std::vector<double> solution_values(n_q_points);
for( const auto& cell : dofHandler.active_cell_iterators() ) {
fe_values.reinit(cell);
fe_values.get_function_values(old_solution, solution_values);
for( unsigned int i=0; i<n_q_points; ++i ) {
*// these values are non-zero, as expected*
* std::cout << "self loop value: " << solution_values[i] << std::endl;*
}
}
Correctly prints the values that I expect.
On Wednesday, January 29, 2020 at 3:14:10 PM UTC-5, Andrew Davis wrote:
>
> I'm trying to implement a time-dependent solver that assembles the system
> using the MeshWorker::loop tool. However, I cannot figure out how get the
> values of the solution at the previous timestep at each quadrature point.
>
> Currently, I store the old solution in a class called "AdvectionProblem"
> such that---I've left out some code in the interest of brevity.
>
> template<unsigned int dim>
> class AdvectionProblem {
> public:
> ...
> void AssembleSystem()
> private:
> ...
> Vector<double> old_solution;
> DoFHandler<dim> dofHandler;
> };
>
> Since we need the old_solution values to assemble the system, I create a
> CellIntegrator class that MeshWorker::loop can use to integrate over each
> cell (I also create similar objects for faces/boundaries):
>
> template<unsigned int dim>
> class IntegrateCell {
> public:
> IntegrateCell(dealii::Vector<double> const& old_solution);
>
> virtual ~IntegrateCell() = default;
>
> void operator()(dealii::MeshWorker::DoFInfo<dim>& dinfo,
> dealii::MeshWorker::IntegrationInfo<dim>& info) const;
>
> private:
>
> /// Store a constant reference to the old solution
> const dealii::Vector<double>& old_solution;
> };
>
> template<unsigned int dim>
> IntegrateCell<dim>::IntegrateCell(Vector<double> const& old_solution) :
> old_solution(old_solution) {}
>
> template<unsigned int dim>
> void IntegrateCell<dim>::operator()(MeshWorker::DoFInfo<dim>& dinfo,
> MeshWorker::IntegrationInfo<dim>& info) const {
> // get some useful objects
> const FEValuesBase<dim>& feValues = info.fe_values();
> FullMatrix<double>& localMatrix = dinfo.matrix(0).matrix;
> Vector<double>& localVector = dinfo.vector(0).block(0);
> const std::vector<double>& JxW = feValues.get_JxW_values();
>
> // get the values of the old solution
> std::vector<double> old_solution_values(feValues.n_quadrature_points);
> feValues.get_function_values(old_solution, old_solution_values);
>
>
> *// these are non-zero, IngtegrateCell is correctly getting the
> old_solution std::cout << old_solution << std::endl;*
>
> // loop over the quadrature points and compute the advection vector at
> the current point
> for( unsigned int point=0; point<feValues.n_quadrature_points; ++point )
> {
>
> *// these are always all zero*
> * std::cout << "value: " << old_solution_values[point] << std::endl;*
> }
> }
>
> I then try to assemble the system using the MeshWorker::loop tool:
>
> template<unsigned int dim>
> void AdvectionProblem<dim>::AssembleSystem() {
> // an object that knows about data structures and local integration
> MeshWorker::IntegrationInfoBox<dim> infoBox;
>
> // initialize the quadrature formula
> const unsigned int nGauss = dofHandler.get_fe().degree + 1;
> infoBox.initialize_gauss_quadrature(nGauss, nGauss, nGauss);
>
> // these are the values we need to integrate our system
> infoBox.initialize_update_flags();
> const UpdateFlags updateFlags = update_quadrature_points | update_values
> | update_gradients | update_JxW_values;
> infoBox.add_update_flags(updateFlags, true, true, true, true);
>
> // initialize the finite element values
> infoBox.initialize(fe, mapping);
>
> // create an object that forwards local data to the assembler
> MeshWorker::DoFInfo<dim> dofInfo(dofHandler);
>
> // create teh assembler---tell it where to put local data
> MeshWorker::Assembler::SystemSimple<SparseMatrix<double>, Vector<double>
> > assembler;
> assembler.initialize(systemMatrix, rhs);
>
> IntegrateCell<dim> cellIntegration(old_solution);
>
> *MeshWorker::loop<dim, dim, MeshWorker::DoFInfo<dim>,
> MeshWorker::IntegrationInfoBox<dim> >(dofHandler.begin_active(),
> dofHandler.end(), dofInfo, infoBox, cellIntegration, nullptr, nullptr,
> assembler);*
> }
>
> *For some reason the feValues.get_function_values(old_solution,
> old_solution_values); call in the CellIntegrator function always sets the
> old_solution_values to zero. Does anyone see what I'm doing wrong?*
>
>
> Thanks!
>
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/4ea6059f-177c-4a7c-8462-5a5482b31b38%40googlegroups.com.
