New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Extend DataOut to complex-valued vectors #1894

Open
bangerth opened this Issue Nov 21, 2015 · 4 comments

Comments

2 participants
@bangerth
Member

bangerth commented Nov 21, 2015

The problem with DataOut and complex valued solution vectors is that (i) it somehow needs to rip real and imaginary parts apart, given that output formats all only support outputting real data; (ii) that DataOut stores the solution vectors it is supposed to output in opaque types that hide whether they are deal.II, PETSc or Trilinos vectors, and whether their elements are floats, doubles, or complex-valued.

To address this, here are three necessary steps:

  • The opaque data types that store data vectors need to be able to tell whether they store a real or complex-valued vector. To this end, let's add a virtual function is_complex_valued() to class internal::DataOut::DataEntryBase that returns whether a data vector is complex-valued. (Fixed by #5333.)
  • Let's add a boolean argument to the various internal::DataOut::DataEntryBase::get_* functions that indicates whether we are asking for the real or imaginary parts. (Fixed by #5335.)
  • In DataOut::build_one_patch, where we call this, we need to ask whether the vector is complex-valued or not (using the function introduced above) and if so, treat real and imaginary parts separately, using the mechanism outlined in part 2 above. (Fixed by #5419.)
  • The same has to be applied to DataOutFaces, DataOutRotation, and DataOutStack.
  • In all of these classes, one also has to deal with the case of the postprocessor and figure out how to sort components in that case.
  • Add instantiations of DataOut_DoFData for complex-valued vectors. (Fixed by #5398.)
  • Add a changelog entry
@bangerth

This comment has been minimized.

Show comment
Hide comment
@bangerth

bangerth Feb 19, 2016

Member

I should add that all of the work above is necessary because we have to separate real and imaginary parts of complex numbers once we want to create graphical output -- no graphical output file I know of understands complex numbers, so we will have to output them as a pair of real values.

Member

bangerth commented Feb 19, 2016

I should add that all of the work above is necessary because we have to separate real and imaginary parts of complex numbers once we want to create graphical output -- no graphical output file I know of understands complex numbers, so we will have to output them as a pair of real values.

@bangerth

This comment has been minimized.

Show comment
Hide comment
@bangerth

bangerth Jul 20, 2017

Member

Here's a testcase:

// ---------------------------------------------------------------------
//
// Copyright (C) 2003 - 2017 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------


// check DataOut for complex vectors

#include "../tests.h"
#include <deal.II/base/logstream.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/block_vector.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/numerics/data_out.h>

#include <fstream>
#include <iomanip>
#include <string>



template <int dim>
void
check ()
{
  Triangulation<dim> tria;
  GridGenerator::hyper_cube(tria, 0., 1.);
  tria.refine_global (1);
  tria.begin_active()->set_refine_flag();
  tria.execute_coarsening_and_refinement ();

  FE_Q<dim> fe(1);
  DoFHandler<dim> dof_handler (tria);
  dof_handler.distribute_dofs (fe);

  Vector<std::complex<double> > v (dof_handler.n_dofs());
  for (unsigned int i=0; i<v.size(); ++i)
    v(i) = std::complex<double>(i,-i);

  DataOut<dim> data_out;
  data_out.attach_dof_handler (dof_handler);
  data_out.add_data_vector (v, "node_data");
  data_out.build_patches ();

  data_out.write_gnuplot (deallog.get_file_stream());
}



int
main()
{
  try
    {
      std::ofstream logfile("output");
      deallog << std::setprecision (2);
      logfile << std::setprecision (2);
      deallog.attach(logfile);
      deallog.depth_console(0);
      deallog.threshold_double(1.e-10);

      check<1> ();
      check<2> ();
      check<3> ();
    }
  catch (std::exception &exc)
    {
      deallog << std::endl << std::endl
              << "----------------------------------------------------"
              << std::endl;
      deallog << "Exception on processing: " << std::endl
              << exc.what() << std::endl
              << "Aborting!" << std::endl
              << "----------------------------------------------------"
              << std::endl;
      return 1;
    }
  catch (...)
    {
      deallog << std::endl << std::endl
              << "----------------------------------------------------"
              << std::endl;
      deallog << "Unknown exception!" << std::endl
              << "Aborting!" << std::endl
              << "----------------------------------------------------"
              << std::endl;
      return 1;
    }
}
Member

bangerth commented Jul 20, 2017

Here's a testcase:

// ---------------------------------------------------------------------
//
// Copyright (C) 2003 - 2017 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------


// check DataOut for complex vectors

#include "../tests.h"
#include <deal.II/base/logstream.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/block_vector.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/numerics/data_out.h>

#include <fstream>
#include <iomanip>
#include <string>



template <int dim>
void
check ()
{
  Triangulation<dim> tria;
  GridGenerator::hyper_cube(tria, 0., 1.);
  tria.refine_global (1);
  tria.begin_active()->set_refine_flag();
  tria.execute_coarsening_and_refinement ();

  FE_Q<dim> fe(1);
  DoFHandler<dim> dof_handler (tria);
  dof_handler.distribute_dofs (fe);

  Vector<std::complex<double> > v (dof_handler.n_dofs());
  for (unsigned int i=0; i<v.size(); ++i)
    v(i) = std::complex<double>(i,-i);

  DataOut<dim> data_out;
  data_out.attach_dof_handler (dof_handler);
  data_out.add_data_vector (v, "node_data");
  data_out.build_patches ();

  data_out.write_gnuplot (deallog.get_file_stream());
}



int
main()
{
  try
    {
      std::ofstream logfile("output");
      deallog << std::setprecision (2);
      logfile << std::setprecision (2);
      deallog.attach(logfile);
      deallog.depth_console(0);
      deallog.threshold_double(1.e-10);

      check<1> ();
      check<2> ();
      check<3> ();
    }
  catch (std::exception &exc)
    {
      deallog << std::endl << std::endl
              << "----------------------------------------------------"
              << std::endl;
      deallog << "Exception on processing: " << std::endl
              << exc.what() << std::endl
              << "Aborting!" << std::endl
              << "----------------------------------------------------"
              << std::endl;
      return 1;
    }
  catch (...)
    {
      deallog << std::endl << std::endl
              << "----------------------------------------------------"
              << std::endl;
      deallog << "Unknown exception!" << std::endl
              << "Aborting!" << std::endl
              << "----------------------------------------------------"
              << std::endl;
      return 1;
    }
}
@bangerth

This comment has been minimized.

Show comment
Hide comment
@bangerth

bangerth Nov 1, 2017

Member

Related in some sene to #2033.

Member

bangerth commented Nov 1, 2017

Related in some sene to #2033.

@bangerth

This comment has been minimized.

Show comment
Hide comment
@bangerth

bangerth Nov 8, 2017

Member

The testcase above now works, once #5419 is merged.

Member

bangerth commented Nov 8, 2017

The testcase above now works, once #5419 is merged.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment