Dear Wolfgang,
I do apologise for that mistake. I hope that what I attach now is OK.
Many thanks once again for any help.
Best wishes,
Katie Leonard
DPhil student in Computational Biology,
The University of Oxford.
________________________________________
From: Wolfgang Bangerth [[email protected]]
Sent: 10 November 2011 17:03
To: Katie Leonard
Cc: Guido Kanschat; [email protected]
Subject: Re: [deal.II] Projecting Initial Values without interpolation
Katie,
> I attach a minimal working program (sorry its still quite lengthy - a lot of
> it is outputting the solution and includes etc.). I have only run this in 2D
> at the moment (its currently set to run in 2D), so I'm not sure that it is
> working properly in 3D - just in case!
That's not what I mean with a minimal program. It takes a long time to
understand someone's code. "Minimal" means: throw out everything you
don't need. For example, if it's the projection of initial values, throw
out everything that has to do with assembling and solving the linear
system -- all I want to see is the way you project initial values and
then output the results to file.
Best
W.
------------------------------------------------------------------------
Wolfgang Bangerth email: [email protected]
www: http://www.math.tamu.edu/~bangerth/
// K H L Leonard modified from step-31.cc
#include <base/quadrature_lib.h>
#include <base/logstream.h>
#include <base/utilities.h>
#include <lac/constraint_matrix.h>
#include <grid/tria.h>
#include <grid/grid_generator.h>
#include <grid/tria_accessor.h>
#include <grid/tria_iterator.h>
#include <grid/tria_boundary_lib.h>
#include <grid/grid_tools.h>
#include <grid/grid_refinement.h>
#include <dofs/dof_handler.h>
#include <dofs/dof_renumbering.h>
#include <dofs/dof_accessor.h>
#include <dofs/dof_tools.h>
#include <fe/fe_q.h>
#include <fe/fe_system.h>
#include <fe/fe_values.h>
#include <numerics/vectors.h>
#include <numerics/data_out.h>
#include <numerics/error_estimator.h>
#include <fstream>
#include <sstream>
#include <iostream>
using namespace dealii;
namespace EquationData
{
template <int dim>
class InitialValues : public Function<dim>
{
public:
InitialValues () : Function<dim>(1) {}
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
virtual void vector_value (const Point<dim> &p,
Vector<double> &value) const;
};
template <int dim>
double
InitialValues<dim>::value (const Point<dim> &p,
const unsigned int component) const
{
if(component == 0)
{
if(p[1]>=-0.2 && p[1] <=0.2 && p[0] <= 0.2-0.5*0.5)
{
return 0.03;
}
else
{
return 0;
}
}
else
{
return 0;
}
}
template <int dim>
void
InitialValues<dim>::vector_value (const Point<dim> &p,
Vector<double> &values) const
{
for (unsigned int c=0; c<this->n_components; ++c)
values(c) = InitialValues<dim>::value (p, c);
}
}
template <int dim>
class ExampleProblem
{
public:
ExampleProblem ();
void run ();
private:
void setup_dofs ();
void output_results () const;
Triangulation<dim> triangulation;
const unsigned int degree;
FESystem<dim> fe;
DoFHandler<dim> dof_handler;
ConstraintMatrix contraints;
Vector<double> solution;
};
template <int dim>
ExampleProblem<dim>::ExampleProblem ()
:
triangulation (Triangulation<dim>::maximum_smoothing),
degree (2),
fe (FE_Q<dim>(degree), 1),
dof_handler (triangulation)
{}
template <int dim>
void ExampleProblem<dim>::setup_dofs ()
{
dof_handler.distribute_dofs (fe);
solution.reinit (dof_handler.n_dofs());
contraints.close();
}
template <int dim>
void ExampleProblem<dim>::output_results () const
{
DataOut<dim> data_out;
data_out.attach_dof_handler (dof_handler);
data_out.add_data_vector (solution, "solution");
data_out.build_patches ();
std::ostringstream filename;
filename << "Example_problem"
<< ".gpl";
std::ofstream output (filename.str().c_str());
data_out.write_gnuplot(output);
}
template <int dim>
void ExampleProblem<dim>::run ()
{
const unsigned int initial_refinement = (dim == 2 ? 7: 2);
GridGenerator::cylinder(triangulation, 0.5, 2);
triangulation.refine_global (initial_refinement);
setup_dofs();
VectorTools::project (dof_handler,
contraints,
QGauss<dim>(degree+2),
EquationData::InitialValues<dim>(),
solution);
output_results ();
}
int main (int argc, char *argv[])
{
try
{
deallog.depth_console (0);
Utilities::System::MPI_InitFinalize mpi_initialization (argc, argv);
ExampleProblem<2> flow_problem;
flow_problem.run ();
}
catch (std::exception &exc)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Exception on processing: " << std::endl
<< exc.what() << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
catch (...)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Unknown exception!" << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
return 0;
}
_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
