Hello all,
I am writing a mesh to a file in ucd format using the GridOut API as follows:
GridOut grid_out;
GridOutFlags::Ucd ucd_flags;
ucd_flags.write_faces = true;
ucd_flags.write_lines = true;
grid_out.set_flags(ucd_flags);
...
std::ofstream grid_ucd_outstream("./grids/optimization_grid_1D.ucd");
std::ofstream grid_gnuplot_outstream("./grids/optimization_grid_1D.gnuplot");
grid_out.write_ucd (triangulation, grid_ucd_outstream);
In another program, I read in the grid like this:
GridIn<dim> grid_in;
if (dim == 2)
{
grid_in.attach_triangulation (triangulation);
std::ifstream input_grid_file("./grids/optimization_grid_2D.ucd");
grid_in.read_ucd(input_grid_file);
}
else if (dim == 3)
{
grid_in.attach_triangulation (triangulation);
std::ifstream input_grid_file("./grids/optimization_grid_3D.ucd");
grid_in.read_ucd(input_grid_file);
}
else
assert (false);
//then I setup the DoFHandler and FE and all other objects, do some time steps
in my solver.
However, the solver does not function correctly anymore. After a single time
step, I get oscillations in the numerical solution. These oscillations do not
occur on an identical grid that is generated in the code by several local
refinement operations (this is my basic sanity check). Moreover, when I try to
refine the grid using DerivativeApproximation::approximate_gradient, i get the
following runtime error:
terminate called after throwing an instance of
'dealii::DerivativeApproximation::ExcInsufficientDirections'
what(): --------------------------------------------------------
An error occurred in line <859> of file
<source/numerics/derivative_approximation.cc> in function
static void dealii::DerivativeApproximation::approximate_cell(const
dealii::Mapping<dim, spacedim>&, const DH<dim, spacedim>&, const InputVector&,
unsigned int, const typename DH<dim, spacedim>::active_cell_iterator&, typename
DerivativeDescription::Derivative&) [with DerivativeDescription =
dealii::DerivativeApproximation::Gradient<2>, int dim = 2, DH =
dealii::DoFHandler, InputVector = dealii::Vector<double>, int spacedim = 2]
The violated condition was:
determinant(Y) != 0
The name and call sequence of the exception was:
ExcInsufficientDirections()
Additional Information:
(none)
Any idea why this happens? Again, if I use a grid that I get from
GridGenerator:hyper_cube (after some local refinements and coarsenings) and is
"identical" to the one I read from the file, everything works OK.
Thanks! and have a good weekend....
-- Mihai
_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
