[cig-commits] commit by bangerth to /var/svn/dealii/aspect
dealii.demon at gmail.com
dealii.demon at gmail.com
Tue Aug 13 08:22:48 PDT 2013
Revision 1837
Split output.h into a file with just the base class and a .cc file that handles and encapsulates the implementation of the various writers.
U trunk/aspect/include/aspect/particle/output.h
A trunk/aspect/source/particle/
A trunk/aspect/source/particle/output.cc
U trunk/aspect/source/postprocess/tracer.cc
http://www.dealii.org/websvn/revision.php?repname=Aspect+Repository&path=%2F&rev=1837&peg=1837
Diff:
Modified: trunk/aspect/include/aspect/particle/output.h
===================================================================
--- trunk/aspect/include/aspect/particle/output.h 2013-08-13 15:16:24 UTC (rev 1836)
+++ trunk/aspect/include/aspect/particle/output.h 2013-08-13 15:22:12 UTC (rev 1837)
@@ -22,12 +22,9 @@
#ifndef __aspect__particle_output_h
#define __aspect__particle_output_h
-#include <deal.II/numerics/data_out.h>
+#include <deal.II/base/mpi.h>
#include <aspect/particle/particle.h>
-#ifdef DEAL_II_HAVE_HDF5
-# include <hdf5.h>
-#endif
namespace aspect
{
@@ -64,8 +61,6 @@
file_index (0)
{}
- virtual ~Output() {};
-
unsigned int self_rank()
{
return Utilities::MPI::this_mpi_process(communicator);
@@ -90,421 +85,21 @@
const double ¤t_time) = 0;
};
- // Blank class to avoid output of data, for example, if particles are used to represent
- // other physical phenomenon that influences the simulation and we don't care about their positions
- template <int dim, class T>
- class NullOutput : public Output<dim, T>
- {
- public:
- virtual std::string output_particle_data(const std::multimap<LevelInd, T> &particles, const double ¤t_time)
- {
- return "";
- };
- };
+ /**
+ * Create an output object given the specified name.
+ */
template <int dim, class T>
- class ASCIIOutput : public Output<dim, T>
- {
- public:
- virtual std::string output_particle_data(const std::multimap<LevelInd, T> &particles,
- const double ¤t_time)
- {
- typename std::multimap<LevelInd, T>::const_iterator it;
- unsigned int i;
- std::string output_file_prefix, output_path_prefix, full_filename;
- std::vector<MPIDataInfo> data_info;
- std::vector<MPIDataInfo>::iterator dit;
- char *particle_data, *p;
+ Output<dim, T> *
+ create_output_object (const std::string &data_format_name);
- output_file_prefix = "particle-" + Utilities::int_to_string (this->file_index, 5);
- output_path_prefix = this->output_dir + output_file_prefix;
- full_filename = output_path_prefix + "." + Utilities::int_to_string(this->self_rank(), 4) + ".txt";
- std::ofstream output (full_filename.c_str());
- if (!output)
- std::cout << "ERROR: proc " << this->self_rank() << " could not create " << full_filename << std::endl;
- // Get the data types
- T::add_mpi_types(data_info);
-
- // Print the header line
- output << "# ";
- for (dit=data_info.begin(); dit!=data_info.end(); ++dit)
- {
- // If it's a 1D element, print just the name, otherwise use []
- if (dit->_num_elems == 1)
- {
- output << dit->_name << " ";
- }
- else
- {
- for (i=0; i<dit->_num_elems; ++i) output << dit->_name << "[" << i << "] ";
- }
- }
- output << "
";
-
- // And print the data for each particle
- particle_data = new char[T::data_len(HDF5_DATA)];
- for (it=particles.begin(); it!=particles.end(); ++it)
- {
- it->second.write_data(HDF5_DATA, particle_data);
- p = particle_data;
- for (dit=data_info.begin(); dit!=data_info.end(); ++dit)
- {
- // Currently assumes all data is double, may need to change this
- for (i=0; i<dit->_num_elems; ++i)
- {
- output << ((double *)p)[0] << " ";
- p += sizeof(double);
- }
- }
- output << "
";
- }
- delete[] particle_data;
-
- output.close();
-
- this->file_index++;
-
- return output_path_prefix;
- };
- };
-
- template <int dim, class T>
- class VTUOutput : public Output<dim, T>
- {
- /**
- * A list of pairs (time, pvtu_filename) that have so far been written
- * and that we will pass to DataOutInterface::write_pvd_record
- * to create a master file that can make the association
- * between simulation time and corresponding file name (this
- * is done because there is no way to store the simulation
- * time inside the .pvtu or .vtu files).
- */
-//TODO: This needs to be serialized
- std::vector<std::pair<double,std::string> > times_and_pvtu_file_names;
- std::vector<std::string> vtu_file_names;
-
- public:
- virtual std::string output_particle_data(const std::multimap<LevelInd, T> &particles, const double ¤t_time)
- {
- typename std::multimap<LevelInd, T>::const_iterator it;
- unsigned int d, i, num_particles;
- DataOut<dim> data_out;
- std::vector<MPIDataInfo> data_info;
- std::vector<MPIDataInfo>::iterator dit;
- unsigned int data_offset;
- char *particle_data, *p;
-
- const std::string output_file_prefix = "particles-" + Utilities::int_to_string (this->file_index, 5);
- const std::string output_path_prefix = this->output_dir + output_file_prefix;
-
- const std::string filename = (output_file_prefix +
- "." +
- Utilities::int_to_string(this->self_rank(), 4) +
- ".vtu");
- const std::string full_filename = (this->output_dir + filename);
-
- std::ofstream output (full_filename.c_str());
- if (!output)
- std::cout << "ERROR: proc "
- << this->self_rank() << " could not create " << filename << std::endl;
-
- num_particles = particles.size();
-
- // Write VTU file XML
- output << "<?xml version=\"1.0\"?>
";
- output << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">
";
- output << " <UnstructuredGrid>
";
- output << " <Piece NumberOfPoints=\"" << num_particles << "\" NumberOfCells=\"" << num_particles << "\">
";
-
- // Go through the particles on this domain and print the position of each one
- output << " <Points>
";
- output << " <DataArray name=\"Position\" type=\"Float64\" NumberOfComponents=\"3\" Format=\"ascii\">
";
- for (it=particles.begin(); it!=particles.end(); ++it)
- {
- output << " " << it->second.location();
-
- // pad with zeros since VTU format wants x/y/z coordinates
- for (d=dim; d<3; ++d)
- output << " 0.0";
-
- output << "
";
- }
- output << " </DataArray>
";
- output << " </Points>
";
-
- // Write cell related data (empty)
- output << " <Cells>
";
- output << " <DataArray type=\"Int32\" Name=\"connectivity\" Format=\"ascii\">
";
- for (i=1; i<=num_particles; ++i)
- output << " " << i-1 << "
";
- output << " </DataArray>
";
- output << " <DataArray type=\"Int32\" Name=\"offsets\" Format=\"ascii\">
";
- for (i=1; i<=num_particles; ++i)
- output << " " << i << "
";
- output << " </DataArray>
";
- output << " <DataArray type=\"UInt8\" Name=\"types\" Format=\"ascii\">
";
- for (i=1; i<=num_particles; ++i)
- output << " 1
";
- output << " </DataArray>
";
- output << " </Cells>
";
-
- // Get the data types
- T::add_mpi_types(data_info);
-
- // Write data for each particle (id, velocity, etc)
- output << " <PointData Scalars=\"scalars\">
";
-
- // Print the data associated with the particles, skipping the first entry (assumed to be position)
- particle_data = new char[T::data_len(HDF5_DATA)];
- dit = data_info.begin();
- data_offset = dit->_num_elems*dit->_elem_size_bytes;
- dit++;
- for (; dit!=data_info.end(); ++dit)
- {
- output << " <DataArray type=\"Float64\" Name=\"" << dit->_name << "\" NumberOfComponents=\"" << (dit->_num_elems == 2 ? 3 : dit->_num_elems) << "\" Format=\"ascii\">
";
- for (it=particles.begin(); it!=particles.end(); ++it)
- {
- it->second.write_data(HDF5_DATA, particle_data);
- p = particle_data+data_offset;
- output << " ";
- for (d=0; d<dit->_num_elems; ++d)
- {
- output << ((double *)p)[0] << " ";
- p += sizeof(double);
- }
- if (dit->_num_elems == 2)
- output << "0 ";
- output << "
";
- }
- data_offset += dit->_num_elems*dit->_elem_size_bytes;
- output << " </DataArray>
";
- }
- delete[] particle_data;
- output << " </PointData>
";
-
- output << " </Piece>
";
- output << " </UnstructuredGrid>
";
- output << "</VTKFile>
";
-
- output.close();
-
-
- // Write the parallel pvtu and pvd files on the root process
- if (this->self_rank() == 0)
- {
- const std::string pvtu_filename = (output_path_prefix + ".pvtu");
-
- std::ofstream pvtu_output (pvtu_filename.c_str());
- if (!pvtu_output) std::cout << "ERROR: Could not create " << filename << std::endl;
-
- pvtu_output << "<?xml version=\"1.0\"?>
";
- pvtu_output << "<VTKFile type=\"PUnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">
";
- pvtu_output << " <PUnstructuredGrid GhostLevel=\"0\">
";
- pvtu_output << " <PPoints>
";
- pvtu_output << " <PDataArray type=\"Float64\" NumberOfComponents=\"3\" format=\"ascii\"/>
";
- pvtu_output << " </PPoints>
";
- pvtu_output << " <PPointData Scalars=\"scalars\">
";
-
- for (dit=data_info.begin()+1; dit!=data_info.end(); ++dit)
- {
- pvtu_output << " <PDataArray type=\"Float64\" Name=\"" << dit->_name << "\" NumberOfComponents=\"" << (dit->_num_elems == 2 ? 3 : dit->_num_elems) << "\" format=\"ascii\"/>
";
- }
- pvtu_output << " </PPointData>
";
- for (i=0; i<this->world_size(); ++i)
- {
- pvtu_output << " <Piece Source=\"" << output_file_prefix << "." << Utilities::int_to_string(i, 4) << ".vtu\"/>
";
- }
- pvtu_output << " </PUnstructuredGrid>
";
- pvtu_output << "</VTKFile>
";
- pvtu_output.close();
-
- times_and_pvtu_file_names.push_back(std::make_pair(current_time,
- pvtu_filename));
- vtu_file_names.push_back (full_filename);
-
- // write .pvd and .visit records
- const std::string pvd_master_filename = (this->output_dir + "particles.pvd");
- std::ofstream pvd_master (pvd_master_filename.c_str());
- data_out.write_pvd_record (pvd_master, times_and_pvtu_file_names);
-
-//TODO: write a global .visit record. this needs a variant of the write_visit_record
-// function
-/*
- const std::string visit_master_filename = (this->output_dir + "particles.visit");
- std::ofstream visit_master (visit_master_filename.c_str());
- data_out.write_visit_record (visit_master, vtu_file_names);
-*/
- }
- this->file_index++;
-
- return output_path_prefix;
- };
- };
-
- template <int dim, class T>
- class HDF5Output : public Output<dim, T>
- {
- // A set of XDMF data objects to create the XDMF file for particles
- std::vector<XDMFEntry> xdmf_entries;
-
- public:
- virtual std::string output_particle_data(const std::multimap<LevelInd, T> &particles, const double ¤t_time)
- {
- typename std::multimap<LevelInd, T>::const_iterator it;
- std::string output_file_prefix, output_path_prefix, full_filename;
-
- output_file_prefix = "particle-" + Utilities::int_to_string (this->file_index, 5);
- output_path_prefix = this->output_dir + output_file_prefix;
-#ifdef DEAL_II_HAVE_HDF5
- // TODO: error checking for H5 calls
- hid_t h5_file_id;
- hid_t plist_id, xfer_plist_id, position_dataset_id, velocity_dataset_id, pid_dataset_id;
- hid_t file_dataspace_id, pos_file_dataspace_id, vel_file_dataspace_id, pid_file_dataspace_id;
- hid_t dim_dataspace_id, one_dim_ds_id;
- hid_t dim_mem_ds_id, one_dim_mem_ds_id;
- hid_t pattr_id;
- herr_t status;
- hsize_t dims[2], offset[2], count[2];
- unsigned int mpi_offset, mpi_count, local_particle_count, global_particle_count, d, i;
-
- double *pos_data, *vel_data, *id_data;
-
- std::string h5_filename = output_path_prefix+".h5";
-
- // Create parallel file access
- plist_id = H5Pcreate(H5P_FILE_ACCESS);
-#ifdef H5_HAVE_PARALLEL
- H5Pset_fapl_mpio(plist_id, this->communicator, MPI_INFO_NULL);
-#endif
-
- // Create the file
- h5_file_id = H5Fcreate(h5_filename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, plist_id);
-
- // Create the file dataspace descriptions
- local_particle_count = particles.size();
- // TODO: error checking on MPI call
- MPI_Comm com = this->communicator;
- MPI_Allreduce(&local_particle_count, &global_particle_count, 1, MPI_UNSIGNED, MPI_SUM, com);
- dims[0] = global_particle_count;
- dims[1] = 3;
- one_dim_ds_id = H5Screate_simple(1, dims, NULL);
- dim_dataspace_id = H5Screate_simple(2, dims, NULL);
-
- // Create the datasets
-#if H5Dcreate_vers == 1
- position_dataset_id = H5Dcreate(h5_file_id, "nodes", H5T_NATIVE_DOUBLE, dim_dataspace_id, H5P_DEFAULT);
- velocity_dataset_id = H5Dcreate(h5_file_id, "velocity", H5T_NATIVE_DOUBLE, dim_dataspace_id, H5P_DEFAULT);
- pid_dataset_id = H5Dcreate(h5_file_id, "id", H5T_NATIVE_DOUBLE, one_dim_ds_id, H5P_DEFAULT);
-#else
- position_dataset_id = H5Dcreate(h5_file_id, "nodes", H5T_NATIVE_DOUBLE, dim_dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
- velocity_dataset_id = H5Dcreate(h5_file_id, "velocity", H5T_NATIVE_DOUBLE, dim_dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
- pid_dataset_id = H5Dcreate(h5_file_id, "id", H5T_NATIVE_DOUBLE, one_dim_ds_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
-#endif
-
- // Close the file dataspaces
- H5Sclose(dim_dataspace_id);
- H5Sclose(one_dim_ds_id);
-
- // Just in case we forget what they are
- count[0] = 1;
- file_dataspace_id = H5Screate_simple (1, count, NULL);
-#if H5Acreate_vers == 1
- pattr_id = H5Acreate(h5_file_id, "Ermahgerd! Pertecrs!", H5T_NATIVE_INT, file_dataspace_id, H5P_DEFAULT);
-#else
- pattr_id = H5Acreate(h5_file_id, "Ermahgerd! Pertecrs!", H5T_NATIVE_INT, file_dataspace_id, H5P_DEFAULT, H5P_DEFAULT);
-#endif
- H5Aclose(pattr_id);
- H5Sclose(file_dataspace_id);
-
- // Count the number of particles on this process and get the offset among all processes
- mpi_count = particles.size();
- MPI_Scan(&mpi_count, &mpi_offset, 1, MPI_UNSIGNED, MPI_SUM, this->communicator);
- count[0] = mpi_count;
- count[1] = 3;
- offset[0] = mpi_offset-mpi_count;
- offset[1] = 0;
-
- // Select the appropriate dataspace for this process
- one_dim_mem_ds_id = H5Screate_simple(1, count, NULL);
- dim_mem_ds_id = H5Screate_simple(2, count, NULL);
- pos_file_dataspace_id = H5Dget_space(position_dataset_id);
- vel_file_dataspace_id = H5Dget_space(velocity_dataset_id);
- pid_file_dataspace_id = H5Dget_space(pid_dataset_id);
-
- // And select the hyperslabs from each dataspace
- status = H5Sselect_hyperslab(pos_file_dataspace_id, H5S_SELECT_SET, offset, NULL, count, NULL);
- status = H5Sselect_hyperslab(vel_file_dataspace_id, H5S_SELECT_SET, offset, NULL, count, NULL);
- status = H5Sselect_hyperslab(pid_file_dataspace_id, H5S_SELECT_SET, offset, NULL, count, NULL);
-
- // Create property list for collective dataset write
- xfer_plist_id = H5Pcreate(H5P_DATASET_XFER);
-#ifdef H5_HAVE_PARALLEL
- H5Pset_dxpl_mpio(xfer_plist_id, H5FD_MPIO_COLLECTIVE);
-#endif
-
- // Read the local particle data
- pos_data = new double[3*particles.size()];
- vel_data = new double[3*particles.size()];
- id_data = new double[particles.size()];
-
- for (i=0,it=particles.begin(); it!=particles.end(); ++i,++it)
- {
- for (d=0; d<dim; ++d)
- {
- pos_data[i*3+d] = it->second.location()(d);
- vel_data[i*3+d] = it->second.velocity()(d);
- }
- if (dim < 3)
- {
- pos_data[i*3+2] = 0;
- vel_data[i*3+2] = 0;
- }
- id_data[i] = it->second.id_num();
- }
-
- // Write particle data to the HDF5 file
- H5Dwrite(position_dataset_id, H5T_NATIVE_DOUBLE, dim_mem_ds_id, pos_file_dataspace_id, xfer_plist_id, pos_data);
- H5Dwrite(velocity_dataset_id, H5T_NATIVE_DOUBLE, dim_mem_ds_id, vel_file_dataspace_id, xfer_plist_id, vel_data);
- H5Dwrite(pid_dataset_id, H5T_NATIVE_DOUBLE, one_dim_mem_ds_id, pid_file_dataspace_id, xfer_plist_id, id_data);
-
- // Clear allocated resources
- delete[] pos_data;
- delete[] vel_data;
- delete[] id_data;
- status = H5Pclose(xfer_plist_id);
- status = H5Sclose(one_dim_mem_ds_id);
- status = H5Sclose(dim_mem_ds_id);
- status = H5Sclose(pos_file_dataspace_id);
- status = H5Sclose(vel_file_dataspace_id);
- status = H5Sclose(pid_file_dataspace_id);
- status = H5Dclose(position_dataset_id);
- status = H5Dclose(velocity_dataset_id);
- status = H5Dclose(pid_dataset_id);
- status = H5Pclose(plist_id);
- status = H5Fclose(h5_file_id);
-
- // Record and output XDMF info on root process
- if (this->self_rank() == 0)
- {
- std::string local_h5_filename = output_file_prefix+".h5";
- XDMFEntry entry(local_h5_filename, current_time, global_particle_count, 0, 3);
- DataOut<dim> data_out;
- const std::string xdmf_filename = (this->output_dir + "particle.xdmf");
-
- entry.add_attribute("velocity", 3);
- entry.add_attribute("id", 1);
- xdmf_entries.push_back(entry);
-
- data_out.write_xdmf_file(xdmf_entries, xdmf_filename.c_str(), this->communicator);
- }
-#endif
-
- this->file_index++;
-
- return output_path_prefix;
- };
- };
+ /**
+ * Return a list of names (separated by '|') of possible writers of graphical
+ * output formats for particle data.
+ */
+ std::string
+ output_object_names ();
}
}
Copied: trunk/aspect/source/particle/output.cc (from rev 1831, trunk/aspect/include/aspect/particle/output.h)
===================================================================
--- trunk/aspect/source/particle/output.cc (rev 0)
+++ trunk/aspect/source/particle/output.cc 2013-08-13 15:22:12 UTC (rev 1837)
@@ -0,0 +1,494 @@
+/*
+ Copyright (C) 2011, 2012 by the authors of the ASPECT code.
+
+ This file is part of ASPECT.
+
+ ASPECT is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2, or (at your option)
+ any later version.
+
+ ASPECT is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with ASPECT; see the file doc/COPYING. If not see
+ <http://www.gnu.org/licenses/>.
+ */
+/* $Id$ */
+
+#include <aspect/particle/output.h>
+#include <deal.II/numerics/data_out.h>
+#include <aspect/particle/particle.h>
+
+#ifdef DEAL_II_HAVE_HDF5
+# include <hdf5.h>
+#endif
+
+namespace aspect
+{
+ namespace Particle
+ {
+ // Blank class to avoid output of data, for example, if particles are used to represent
+ // other physical phenomenon that influences the simulation and we don't care about their positions
+ template <int dim, class T>
+ class NullOutput : public Output<dim, T>
+ {
+ public:
+ virtual std::string output_particle_data(const std::multimap<LevelInd, T> &particles, const double ¤t_time)
+ {
+ return "";
+ };
+ };
+
+ template <int dim, class T>
+ class ASCIIOutput : public Output<dim, T>
+ {
+ public:
+ virtual std::string output_particle_data(const std::multimap<LevelInd, T> &particles,
+ const double ¤t_time)
+ {
+ typename std::multimap<LevelInd, T>::const_iterator it;
+ unsigned int i;
+ std::string output_file_prefix, output_path_prefix, full_filename;
+ std::vector<MPIDataInfo> data_info;
+ std::vector<MPIDataInfo>::iterator dit;
+ char *particle_data, *p;
+
+ output_file_prefix = "particle-" + Utilities::int_to_string (this->file_index, 5);
+ output_path_prefix = this->output_dir + output_file_prefix;
+ full_filename = output_path_prefix + "." + Utilities::int_to_string(this->self_rank(), 4) + ".txt";
+ std::ofstream output (full_filename.c_str());
+ if (!output)
+ std::cout << "ERROR: proc " << this->self_rank() << " could not create " << full_filename << std::endl;
+
+ // Get the data types
+ T::add_mpi_types(data_info);
+
+ // Print the header line
+ output << "# ";
+ for (dit=data_info.begin(); dit!=data_info.end(); ++dit)
+ {
+ // If it's a 1D element, print just the name, otherwise use []
+ if (dit->_num_elems == 1)
+ {
+ output << dit->_name << " ";
+ }
+ else
+ {
+ for (i=0; i<dit->_num_elems; ++i) output << dit->_name << "[" << i << "] ";
+ }
+ }
+ output << "
";
+
+ // And print the data for each particle
+ particle_data = new char[T::data_len(HDF5_DATA)];
+ for (it=particles.begin(); it!=particles.end(); ++it)
+ {
+ it->second.write_data(HDF5_DATA, particle_data);
+ p = particle_data;
+ for (dit=data_info.begin(); dit!=data_info.end(); ++dit)
+ {
+ // Currently assumes all data is double, may need to change this
+ for (i=0; i<dit->_num_elems; ++i)
+ {
+ output << ((double *)p)[0] << " ";
+ p += sizeof(double);
+ }
+ }
+ output << "
";
+ }
+ delete[] particle_data;
+
+ output.close();
+
+ this->file_index++;
+
+ return output_path_prefix;
+ };
+ };
+
+
+
+ template <int dim, class T>
+ class VTUOutput : public Output<dim, T>
+ {
+ /**
+ * A list of pairs (time, pvtu_filename) that have so far been written
+ * and that we will pass to DataOutInterface::write_pvd_record
+ * to create a master file that can make the association
+ * between simulation time and corresponding file name (this
+ * is done because there is no way to store the simulation
+ * time inside the .pvtu or .vtu files).
+ */
+//TODO: This needs to be serialized
+ std::vector<std::pair<double,std::string> > times_and_pvtu_file_names;
+ std::vector<std::string> vtu_file_names;
+
+ public:
+ virtual std::string output_particle_data(const std::multimap<LevelInd, T> &particles, const double ¤t_time)
+ {
+ typename std::multimap<LevelInd, T>::const_iterator it;
+ unsigned int d, i, num_particles;
+ DataOut<dim> data_out;
+ std::vector<MPIDataInfo> data_info;
+ std::vector<MPIDataInfo>::iterator dit;
+ unsigned int data_offset;
+ char *particle_data, *p;
+
+ const std::string output_file_prefix = "particles-" + Utilities::int_to_string (this->file_index, 5);
+ const std::string output_path_prefix = this->output_dir + output_file_prefix;
+
+ const std::string filename = (output_file_prefix +
+ "." +
+ Utilities::int_to_string(this->self_rank(), 4) +
+ ".vtu");
+ const std::string full_filename = (this->output_dir + filename);
+
+ std::ofstream output (full_filename.c_str());
+ if (!output)
+ std::cout << "ERROR: proc "
+ << this->self_rank() << " could not create " << filename << std::endl;
+
+ num_particles = particles.size();
+
+ // Write VTU file XML
+ output << "<?xml version=\"1.0\"?>
";
+ output << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">
";
+ output << " <UnstructuredGrid>
";
+ output << " <Piece NumberOfPoints=\"" << num_particles << "\" NumberOfCells=\"" << num_particles << "\">
";
+
+ // Go through the particles on this domain and print the position of each one
+ output << " <Points>
";
+ output << " <DataArray name=\"Position\" type=\"Float64\" NumberOfComponents=\"3\" Format=\"ascii\">
";
+ for (it=particles.begin(); it!=particles.end(); ++it)
+ {
+ output << " " << it->second.location();
+
+ // pad with zeros since VTU format wants x/y/z coordinates
+ for (d=dim; d<3; ++d)
+ output << " 0.0";
+
+ output << "
";
+ }
+ output << " </DataArray>
";
+ output << " </Points>
";
+
+ // Write cell related data (empty)
+ output << " <Cells>
";
+ output << " <DataArray type=\"Int32\" Name=\"connectivity\" Format=\"ascii\">
";
+ for (i=1; i<=num_particles; ++i)
+ output << " " << i-1 << "
";
+ output << " </DataArray>
";
+ output << " <DataArray type=\"Int32\" Name=\"offsets\" Format=\"ascii\">
";
+ for (i=1; i<=num_particles; ++i)
+ output << " " << i << "
";
+ output << " </DataArray>
";
+ output << " <DataArray type=\"UInt8\" Name=\"types\" Format=\"ascii\">
";
+ for (i=1; i<=num_particles; ++i)
+ output << " 1
";
+ output << " </DataArray>
";
+ output << " </Cells>
";
+
+ // Get the data types
+ T::add_mpi_types(data_info);
+
+ // Write data for each particle (id, velocity, etc)
+ output << " <PointData Scalars=\"scalars\">
";
+
+ // Print the data associated with the particles, skipping the first entry (assumed to be position)
+ particle_data = new char[T::data_len(HDF5_DATA)];
+ dit = data_info.begin();
+ data_offset = dit->_num_elems*dit->_elem_size_bytes;
+ dit++;
+ for (; dit!=data_info.end(); ++dit)
+ {
+ output << " <DataArray type=\"Float64\" Name=\"" << dit->_name << "\" NumberOfComponents=\"" << (dit->_num_elems == 2 ? 3 : dit->_num_elems) << "\" Format=\"ascii\">
";
+ for (it=particles.begin(); it!=particles.end(); ++it)
+ {
+ it->second.write_data(HDF5_DATA, particle_data);
+ p = particle_data+data_offset;
+ output << " ";
+ for (d=0; d<dit->_num_elems; ++d)
+ {
+ output << ((double *)p)[0] << " ";
+ p += sizeof(double);
+ }
+ if (dit->_num_elems == 2)
+ output << "0 ";
+ output << "
";
+ }
+ data_offset += dit->_num_elems*dit->_elem_size_bytes;
+ output << " </DataArray>
";
+ }
+ delete[] particle_data;
+ output << " </PointData>
";
+
+ output << " </Piece>
";
+ output << " </UnstructuredGrid>
";
+ output << "</VTKFile>
";
+
+ output.close();
+
+
+ // Write the parallel pvtu and pvd files on the root process
+ if (this->self_rank() == 0)
+ {
+ const std::string pvtu_filename = (output_path_prefix + ".pvtu");
+
+ std::ofstream pvtu_output (pvtu_filename.c_str());
+ if (!pvtu_output) std::cout << "ERROR: Could not create " << filename << std::endl;
+
+ pvtu_output << "<?xml version=\"1.0\"?>
";
+ pvtu_output << "<VTKFile type=\"PUnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">
";
+ pvtu_output << " <PUnstructuredGrid GhostLevel=\"0\">
";
+ pvtu_output << " <PPoints>
";
+ pvtu_output << " <PDataArray type=\"Float64\" NumberOfComponents=\"3\" format=\"ascii\"/>
";
+ pvtu_output << " </PPoints>
";
+ pvtu_output << " <PPointData Scalars=\"scalars\">
";
+
+ for (dit=data_info.begin()+1; dit!=data_info.end(); ++dit)
+ {
+ pvtu_output << " <PDataArray type=\"Float64\" Name=\"" << dit->_name << "\" NumberOfComponents=\"" << (dit->_num_elems == 2 ? 3 : dit->_num_elems) << "\" format=\"ascii\"/>
";
+ }
+ pvtu_output << " </PPointData>
";
+ for (i=0; i<this->world_size(); ++i)
+ {
+ pvtu_output << " <Piece Source=\"" << output_file_prefix << "." << Utilities::int_to_string(i, 4) << ".vtu\"/>
";
+ }
+ pvtu_output << " </PUnstructuredGrid>
";
+ pvtu_output << "</VTKFile>
";
+ pvtu_output.close();
+
+ times_and_pvtu_file_names.push_back(std::make_pair(current_time,
+ pvtu_filename));
+ vtu_file_names.push_back (full_filename);
+
+ // write .pvd and .visit records
+ const std::string pvd_master_filename = (this->output_dir + "particles.pvd");
+ std::ofstream pvd_master (pvd_master_filename.c_str());
+ data_out.write_pvd_record (pvd_master, times_and_pvtu_file_names);
+
+//TODO: write a global .visit record. this needs a variant of the write_visit_record
+// function
+/*
+ const std::string visit_master_filename = (this->output_dir + "particles.visit");
+ std::ofstream visit_master (visit_master_filename.c_str());
+ data_out.write_visit_record (visit_master, vtu_file_names);
+*/
+ }
+ this->file_index++;
+
+ return output_path_prefix;
+ };
+ };
+
+
+
+ template <int dim, class T>
+ class HDF5Output : public Output<dim, T>
+ {
+ private:
+ // A set of XDMF data objects to create the XDMF file for particles
+ std::vector<XDMFEntry> xdmf_entries;
+
+ public:
+ virtual std::string output_particle_data(const std::multimap<LevelInd, T> &particles, const double ¤t_time)
+ {
+ typename std::multimap<LevelInd, T>::const_iterator it;
+ std::string output_file_prefix, output_path_prefix, full_filename;
+
+ output_file_prefix = "particle-" + Utilities::int_to_string (this->file_index, 5);
+ output_path_prefix = this->output_dir + output_file_prefix;
+#ifdef DEAL_II_HAVE_HDF5
+ // TODO: error checking for H5 calls
+ hid_t h5_file_id;
+ hid_t plist_id, xfer_plist_id, position_dataset_id, velocity_dataset_id, pid_dataset_id;
+ hid_t file_dataspace_id, pos_file_dataspace_id, vel_file_dataspace_id, pid_file_dataspace_id;
+ hid_t dim_dataspace_id, one_dim_ds_id;
+ hid_t dim_mem_ds_id, one_dim_mem_ds_id;
+ hid_t pattr_id;
+ herr_t status;
+ hsize_t dims[2], offset[2], count[2];
+ unsigned int mpi_offset, mpi_count, local_particle_count, global_particle_count, d, i;
+
+ double *pos_data, *vel_data, *id_data;
+
+ std::string h5_filename = output_path_prefix+".h5";
+
+ // Create parallel file access
+ plist_id = H5Pcreate(H5P_FILE_ACCESS);
+#ifdef H5_HAVE_PARALLEL
+ H5Pset_fapl_mpio(plist_id, this->communicator, MPI_INFO_NULL);
+#endif
+
+ // Create the file
+ h5_file_id = H5Fcreate(h5_filename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, plist_id);
+
+ // Create the file dataspace descriptions
+ local_particle_count = particles.size();
+ // TODO: error checking on MPI call
+ MPI_Comm com = this->communicator;
+ MPI_Allreduce(&local_particle_count, &global_particle_count, 1, MPI_UNSIGNED, MPI_SUM, com);
+ dims[0] = global_particle_count;
+ dims[1] = 3;
+ one_dim_ds_id = H5Screate_simple(1, dims, NULL);
+ dim_dataspace_id = H5Screate_simple(2, dims, NULL);
+
+ // Create the datasets
+#if H5Dcreate_vers == 1
+ position_dataset_id = H5Dcreate(h5_file_id, "nodes", H5T_NATIVE_DOUBLE, dim_dataspace_id, H5P_DEFAULT);
+ velocity_dataset_id = H5Dcreate(h5_file_id, "velocity", H5T_NATIVE_DOUBLE, dim_dataspace_id, H5P_DEFAULT);
+ pid_dataset_id = H5Dcreate(h5_file_id, "id", H5T_NATIVE_DOUBLE, one_dim_ds_id, H5P_DEFAULT);
+#else
+ position_dataset_id = H5Dcreate(h5_file_id, "nodes", H5T_NATIVE_DOUBLE, dim_dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+ velocity_dataset_id = H5Dcreate(h5_file_id, "velocity", H5T_NATIVE_DOUBLE, dim_dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+ pid_dataset_id = H5Dcreate(h5_file_id, "id", H5T_NATIVE_DOUBLE, one_dim_ds_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+#endif
+
+ // Close the file dataspaces
+ H5Sclose(dim_dataspace_id);
+ H5Sclose(one_dim_ds_id);
+
+ // Just in case we forget what they are
+ count[0] = 1;
+ file_dataspace_id = H5Screate_simple (1, count, NULL);
+#if H5Acreate_vers == 1
+ pattr_id = H5Acreate(h5_file_id, "Ermahgerd! Pertecrs!", H5T_NATIVE_INT, file_dataspace_id, H5P_DEFAULT);
+#else
+ pattr_id = H5Acreate(h5_file_id, "Ermahgerd! Pertecrs!", H5T_NATIVE_INT, file_dataspace_id, H5P_DEFAULT, H5P_DEFAULT);
+#endif
+ H5Aclose(pattr_id);
+ H5Sclose(file_dataspace_id);
+
+ // Count the number of particles on this process and get the offset among all processes
+ mpi_count = particles.size();
+ MPI_Scan(&mpi_count, &mpi_offset, 1, MPI_UNSIGNED, MPI_SUM, this->communicator);
+ count[0] = mpi_count;
+ count[1] = 3;
+ offset[0] = mpi_offset-mpi_count;
+ offset[1] = 0;
+
+ // Select the appropriate dataspace for this process
+ one_dim_mem_ds_id = H5Screate_simple(1, count, NULL);
+ dim_mem_ds_id = H5Screate_simple(2, count, NULL);
+ pos_file_dataspace_id = H5Dget_space(position_dataset_id);
+ vel_file_dataspace_id = H5Dget_space(velocity_dataset_id);
+ pid_file_dataspace_id = H5Dget_space(pid_dataset_id);
+
+ // And select the hyperslabs from each dataspace
+ status = H5Sselect_hyperslab(pos_file_dataspace_id, H5S_SELECT_SET, offset, NULL, count, NULL);
+ status = H5Sselect_hyperslab(vel_file_dataspace_id, H5S_SELECT_SET, offset, NULL, count, NULL);
+ status = H5Sselect_hyperslab(pid_file_dataspace_id, H5S_SELECT_SET, offset, NULL, count, NULL);
+
+ // Create property list for collective dataset write
+ xfer_plist_id = H5Pcreate(H5P_DATASET_XFER);
+#ifdef H5_HAVE_PARALLEL
+ H5Pset_dxpl_mpio(xfer_plist_id, H5FD_MPIO_COLLECTIVE);
+#endif
+
+ // Read the local particle data
+ pos_data = new double[3*particles.size()];
+ vel_data = new double[3*particles.size()];
+ id_data = new double[particles.size()];
+
+ for (i=0,it=particles.begin(); it!=particles.end(); ++i,++it)
+ {
+ for (d=0; d<dim; ++d)
+ {
+ pos_data[i*3+d] = it->second.location()(d);
+ vel_data[i*3+d] = it->second.velocity()(d);
+ }
+ if (dim < 3)
+ {
+ pos_data[i*3+2] = 0;
+ vel_data[i*3+2] = 0;
+ }
+ id_data[i] = it->second.id_num();
+ }
+
+ // Write particle data to the HDF5 file
+ H5Dwrite(position_dataset_id, H5T_NATIVE_DOUBLE, dim_mem_ds_id, pos_file_dataspace_id, xfer_plist_id, pos_data);
+ H5Dwrite(velocity_dataset_id, H5T_NATIVE_DOUBLE, dim_mem_ds_id, vel_file_dataspace_id, xfer_plist_id, vel_data);
+ H5Dwrite(pid_dataset_id, H5T_NATIVE_DOUBLE, one_dim_mem_ds_id, pid_file_dataspace_id, xfer_plist_id, id_data);
+
+ // Clear allocated resources
+ delete[] pos_data;
+ delete[] vel_data;
+ delete[] id_data;
+ status = H5Pclose(xfer_plist_id);
+ status = H5Sclose(one_dim_mem_ds_id);
+ status = H5Sclose(dim_mem_ds_id);
+ status = H5Sclose(pos_file_dataspace_id);
+ status = H5Sclose(vel_file_dataspace_id);
+ status = H5Sclose(pid_file_dataspace_id);
+ status = H5Dclose(position_dataset_id);
+ status = H5Dclose(velocity_dataset_id);
+ status = H5Dclose(pid_dataset_id);
+ status = H5Pclose(plist_id);
+ status = H5Fclose(h5_file_id);
+
+ // Record and output XDMF info on root process
+ if (this->self_rank() == 0)
+ {
+ std::string local_h5_filename = output_file_prefix+".h5";
+ XDMFEntry entry(local_h5_filename, current_time, global_particle_count, 0, 3);
+ DataOut<dim> data_out;
+ const std::string xdmf_filename = (this->output_dir + "particle.xdmf");
+
+ entry.add_attribute("velocity", 3);
+ entry.add_attribute("id", 1);
+ xdmf_entries.push_back(entry);
+
+ data_out.write_xdmf_file(xdmf_entries, xdmf_filename.c_str(), this->communicator);
+ }
+#endif
+
+ this->file_index++;
+
+ return output_path_prefix;
+ };
+ };
+
+
+
+ template <int dim, class T>
+ Output<dim,T> *
+ create_output_object (const std::string &data_output_format)
+ {
+ if (data_output_format == "ascii")
+ return new Particle::ASCIIOutput<dim,T>();
+ else if (data_output_format == "vtu")
+ return new Particle::VTUOutput<dim,T>();
+ else if (data_output_format == "hdf5")
+ return new Particle::HDF5Output<dim,T>();
+ else if (data_output_format == "none")
+ return new Particle::NullOutput<dim,T>();
+ else
+ Assert (false, ExcNotImplemented());
+
+ return 0;
+ }
+
+
+ std::string
+ output_object_names ()
+ {
+ return ("none|"
+ "ascii|"
+ "vtu|"
+ "hdf5");
+ }
+
+
+ // explicit instantiations
+ template
+ Output<2,Particle::BaseParticle<2> > *
+ create_output_object (const std::string &data_output_format);
+ template
+ Output<3,Particle::BaseParticle<3> > *
+ create_output_object (const std::string &data_output_format);
+ }
+}
Modified: trunk/aspect/source/postprocess/tracer.cc
===================================================================
--- trunk/aspect/source/postprocess/tracer.cc 2013-08-13 15:16:24 UTC (rev 1836)
+++ trunk/aspect/source/postprocess/tracer.cc 2013-08-13 15:22:12 UTC (rev 1837)
@@ -43,22 +43,11 @@
if (!initialized)
{
-//TODO: Use factory methods here
// Create an output object depending on what the parameters specify
- if (data_output_format == "ascii")
- output = new Particle::ASCIIOutput<dim,Particle::BaseParticle<dim> >();
- else if (data_output_format == "vtu")
- output = new Particle::VTUOutput<dim,Particle::BaseParticle<dim> >();
- else if (data_output_format == "hdf5")
- output = new Particle::HDF5Output<dim,Particle::BaseParticle<dim> >();
- else if (data_output_format == "none")
- output = new Particle::NullOutput<dim,Particle::BaseParticle<dim> >();
- else
- Assert (false, ExcNotImplemented());
-
- // Set the output directory for the particle output to be stored in
+ output = Particle::create_output_object<dim,Particle::BaseParticle<dim> > (data_output_format);
output->set_output_directory(this->get_output_directory());
+ //TODO: Use factory methods here
// Create an integrator object depending on the specified parameter
if (integration_scheme == "euler")
integrator = new Particle::EulerIntegrator<dim, Particle::BaseParticle<dim> >;
@@ -153,11 +142,7 @@
"'Use years in output instead of seconds' parameter is set; "
"seconds otherwise.");
prm.declare_entry("Data output format", "none",
- Patterns::Selection("none|"
- "ascii|"
- "vtu|"
- "hdf5"
- ),
+ Patterns::Selection(Particle::output_object_names()),
More information about the CIG-COMMITS
mailing list