Sure, attached is a *.cpp file that should run. The cout statement on line
87 prints non-zero values for `old_solution` but the cout statement online
92 always prints zero.
On Thursday, January 30, 2020 at 10:49:18 AM UTC-5, Wolfgang Bangerth wrote:
>
> On 1/30/20 8:45 AM, Andrew Davis wrote:
> > I thought the same thing---that 'old_solution' would be a zero
> vector---but
> > when I add the line `std::cout << old_solution << std::endl;` in the
> function
> > IntegrateCell<dim>::operator() (just after the `get_function_values'
> call in
> > the original post) prints a non-zero vector that I expect.
>
> That does not sound credible to me. Can you construct a small example that
> does exactly that? I.e., print out, for example, the l2-norm of
> old_solution
> right before the call to fe_values.get_function_values(), and the l2-norm
> of
> old_solution_values right after.
>
> Best
> W.
>
> --
> ------------------------------------------------------------------------
> Wolfgang Bangerth email: bang...@colostate.edu
> <javascript:>
> www: http://www.math.colostate.edu/~bangerth/
>
>
--
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/4dd7f499-1eeb-41b2-bed4-d715354007d7%40googlegroups.com.
#include <iostream>
#include <fstream>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/solver_richardson.h>
#include <deal.II/lac/precondition_block.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/mapping_q1.h>
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/derivative_approximation.h>
#include <deal.II/meshworker/dof_info.h>
#include <deal.II/meshworker/integration_info.h>
#include <deal.II/meshworker/simple.h>
#include <deal.II/meshworker/loop.h>
using namespace dealii;
/// The initial conditions for the mass
template<unsigned int dim>
class InitialConditionMass : public dealii::Function<dim> {
public:
InitialConditionMass() = default;
virtual ~InitialConditionMass() = default;
virtual double value(dealii::Point<dim> const& p, const unsigned int component=0) const override;
private:
};
template<unsigned int dim>
double InitialConditionMass<dim>::value(dealii::Point<dim> const& p, const unsigned int component) const {
(void)component;
Assert(component==0, ExcIndexRange(component, 0, 1));
return p(0);
}
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 {
const FEValuesBase<dim>& feValues = info.fe_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);
//std::cout << "old_solution: " << old_solution.l2_norm() << std::endl;
std::cout << "old_solution: " << 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 ) {
std::cout << "point: " << feValues.quadrature_point(point) << std::endl;
std::cout << "value: " << old_solution_values[point] << std::endl;
}
}
/// The advection problem class
template<unsigned int dim>
class AdvectionProblem {
public:
/// Default constructor
AdvectionProblem();
virtual ~AdvectionProblem() = default;
void SetupSystem();
void AssembleSystem();
void Run();
private:
/// Set up the initial condition
/**
@param[in] vec Set the intial condition onto this vector
*/
void SetInitialCondition(dealii::Vector<double>& vec);
/// Globally refine the mesh
/**
@param[in] nrefine The number of global refinements
*/
void InitializeGrid(unsigned int const nrefine);
using DoFInfo = dealii::MeshWorker::DoFInfo<dim>;
using CellInfo = dealii::MeshWorker::IntegrationInfo<dim>;
dealii::Triangulation<dim> triangulation;
const dealii::MappingQ1<dim> mapping;
dealii::FE_DGQ<dim> fe;
dealii::DoFHandler<dim> dofHandler;
dealii::AffineConstraints<double> constraints;
dealii::SparsityPattern sparsityPattern;
dealii::SparseMatrix<double> systemMatrix;
/// The solution to the advection equation at the current timestep
dealii::Vector<double> solution;
/// The solution to the advection equation at the previous timestep
dealii::Vector<double> old_solution;
dealii::Vector<double> rhs;
/// The length of each timestep
const double timeStep = 1.0/3.0;
/// The current time
double currTime;
/// The current timestep
unsigned int step = 0;
/// The final time
const double finalTime = 1.0;
};
template<unsigned int dim>
AdvectionProblem<dim>::AdvectionProblem() : mapping(), fe(1), dofHandler(triangulation), currTime(timeStep) {
// initialize a square grid and refine global
InitializeGrid(5);
// set up the system on this initial grid
SetupSystem();
// reset the old solution
old_solution.reinit(dofHandler.n_dofs());
}
template<unsigned int dim>
void AdvectionProblem<dim>::SetupSystem() {
// distribute the dofs
dofHandler.distribute_dofs(fe);
// generate a sparsity pattern
DynamicSparsityPattern dsp(dofHandler.n_dofs());
DoFTools::make_flux_sparsity_pattern(dofHandler, dsp);
sparsityPattern.copy_from(dsp);
// set up the strucutre of the components
systemMatrix.reinit(sparsityPattern);
solution.reinit(dofHandler.n_dofs());
rhs.reinit(dofHandler.n_dofs());
constraints.close();
}
template<unsigned int dim>
void AdvectionProblem<dim>::SetInitialCondition(Vector<double>& vec) {
VectorTools::project(dofHandler, constraints, QGauss<dim>(fe.degree+1), InitialConditionMass<dim>(), vec);
}
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 the 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);
}
template<unsigned int dim>
void AdvectionProblem<dim>::InitializeGrid(unsigned int const nrefine) {
GridGenerator::hyper_cube(triangulation);
triangulation.refine_global(nrefine);
}
template<unsigned int dim>
void AdvectionProblem<dim>::Run() {
// project the initial conditions---on whatever grid is currently setup (constructor initializes something)
SetInitialCondition(old_solution);
// loop through time
for( ; currTime<finalTime+timeStep/2.0; currTime+=timeStep, ++step ) {
std::cout << "time: " << currTime << " of " << finalTime << std::endl;
// set up the sparsity pattern for the system and reinit the solution
SetupSystem();
solution = old_solution;
// assemble the system
AssembleSystem();
}
}
int main(int argc, char **argv) {
AdvectionProblem<2> dgmethod;
dgmethod.Run();
}
