OutputManager.cpp
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Thu, Jul 11, 17:41

OutputManager.cpp
View Options

	#include <string>
	#include <fstream>
	#include <stdio.h>
	#include <sstream>
	#include "OutputManager.h"
	#include "ATC_Error.h"

	//#define EXTENDED_ERROR_CHECKING

	namespace ATC {

	static const int kFieldPrecison = 12;
	static const int kFieldWidth = kFieldPrecison + 6;

	static const int kFileNameSize = 26;
	static string tensor_component_names[9] = {"11","12","13",
	"21","22","23",
	"31","32","33"};
	static string sym_tensor_component_names[6] = {"11","22","33","12","13","23"};
	static string vector_component_names[3] = {"_X","_Y","_Z"};
	static string list_component_names[26] = {"_a","_b","_c","_d","_e","_f","_g","_h","_i","_j","_k","_l","_m","_n","_o","_p","_q","_r","_s","_t","_u","_v","_w","_x","_y","_z"};

	string* get_component_names(int type) {
	string* componentNames = list_component_names;
	if (type==VECTOR_OUTPUT)
	componentNames = vector_component_names;
	else if (type == SYM_TENSOR_OUTPUT)
	componentNames = sym_tensor_component_names;
	else if (type == TENSOR_OUTPUT)
	componentNames = tensor_component_names;
	return componentNames;
	}


	//-----------------------------------------------------------------------------
	//*
	//-----------------------------------------------------------------------------
	OutputManager::OutputManager(string outputPrefix, OutputType otype)
	: outputPrefix_(outputPrefix),
	dataType_(POINT),
	outputType_(otype),
	firstStep_(true),
	writeGlobalsHeader_(true),
	firstGlobalsWrite_(true),
	tensorToComponents_(false), // paraview does not support tensors
	vectorToComponents_(false),
	initialized_(false)
	{}
	//-----------------------------------------------------------------------------
	//*
	//-----------------------------------------------------------------------------
	OutputManager::OutputManager()
	: outputPrefix_("NULL"),
	dataType_(POINT),
	outputType_(ENSIGHT),
	firstStep_(true),
	writeGlobalsHeader_(true),
	firstGlobalsWrite_(true),
	tensorToComponents_(false), // paraview does not support tensors
	vectorToComponents_(false),
	initialized_(false)
	{}
	//-----------------------------------------------------------------------------
	//*
	//-----------------------------------------------------------------------------
	OutputManager::~OutputManager() {}
	//-----------------------------------------------------------------------------
	//*
	//-----------------------------------------------------------------------------
	void OutputManager::set_option(OutputOption option, bool value) {
	if (option == OUTPUT_VECTOR_COMPONENTS) vectorToComponents_ = value;
	else if (option == OUTPUT_TENSOR_COMPONENTS) tensorToComponents_ = value;
	else throw ATC_Error(0,"unsupported output option");
	};

	//-----------------------------------------------------------------------------
	//*
	//-----------------------------------------------------------------------------
	void OutputManager::initialize(string outputPrefix, OutputType otype)
	{
	if (outputPrefix_ != outputPrefix ) { // new stream with existing object
	outputPrefix_ = outputPrefix;
	initialized_ = false;
	}
	outputTimes_.clear();
	outputType_ = otype;
	firstStep_ = true;
	firstGlobalsWrite_ = true;
	writeGlobalsHeader_ = true;
	}
	//-----------------------------------------------------------------------------
	//*
	//-----------------------------------------------------------------------------
	// Dump text-based fields to disk for later restart
	void OutputManager::write_restart_file(string fileName, OUTPUT_LIST *data)
	{
	FILE * fp=NULL;
	fp=fopen(fileName.c_str(),"wb"); // open
	OUTPUT_LIST::iterator iter;
	for (iter = data->begin(); iter != data->end(); iter++) {
	const MATRIX* field_data = iter->second;
	for (int i = 0; i < field_data->nRows(); ++i) {
	for (int j = 0; j < field_data->nCols(); ++j) {
	double x = (*field_data)(i,j);
	fwrite(&x,sizeof(double),1,fp);
	}
	}
	}
	fclose(fp);
	}
	//-----------------------------------------------------------------------------
	//*
	//-----------------------------------------------------------------------------
	// Read a file corresponding to a write by write_restart_file
	void OutputManager::read_restart_file(string fileName, OUTPUT_LIST *data)
	{
	FILE * fp=NULL;
	fp=fopen(fileName.c_str(),"rb"); // open
	OUTPUT_LIST::iterator iter;
	for (iter = data->begin(); iter != data->end(); iter++) {
	MATRIX* field_data = iter->second;
	for (int i = 0; i < field_data->nRows(); ++i) {
	for (int j = 0; j < field_data->nCols(); ++j) {
	double myVal;
	fread(&myVal,sizeof(double),1,fp);
	(*field_data)(i,j) = myVal;
	}
	}
	}
	fclose(fp);
	}
	//-----------------------------------------------------------------------------
	//*
	//-----------------------------------------------------------------------------
	void OutputManager::write_globals(void)
	{
	if ( outputPrefix_ == "NULL") return;
	string file = outputPrefix_ + ".GLOBALS";
	ofstream text;
	if ( firstGlobalsWrite_ ) text.open(file.c_str(),ios_base::out);
	else text.open(file.c_str(),ios_base::app);
	firstGlobalsWrite_ = false;

	map<string, double>::iterator iter;
	// header
	if ( firstStep_ \|\| writeGlobalsHeader_) {
	text << "# time:1 ";
	int index = 2;
	for (iter = globalData_.begin(); iter != globalData_.end(); iter++)
	{
	string name = iter->first;
	string str; stringstream out; out << ":" << index++;
	str = out.str();
	name.append(str);
	text.width(kFieldWidth); text << name << " ";
	}
	text << '\n';
	}
	writeGlobalsHeader_ = false;
	// data
	text.width(kFieldWidth); text << outputTimes_[outputTimes_.size()-1] << " ";
	for (iter = globalData_.begin();
	iter != globalData_.end(); iter++) {
	double value = iter->second;
	text << setw(kFieldWidth) << std::scientific << std::setprecision(kFieldPrecison) << value << " ";
	}
	text << "\n";
	}
	//-----------------------------------------------------------------------------
	//*
	//-----------------------------------------------------------------------------
	void OutputManager::write_geometry(const MATRIX &coordinates,
	const Array2D<int> *connectivities)
	{
	if ( outputPrefix_ == "NULL") return;
	// geometry based on a reference configuration
	string geom_file_name = outputPrefix_ + ".geo";
	string geom_file_text = outputPrefix_ + ".XYZ";

	// open file
	FILE * fp=NULL;
	ofstream text;
	char buffer[80];
	if ( ! initialized_ ) {
	fp=fopen(geom_file_name.c_str(),"wb"); // open
	strcpy(buffer,"C Binary");
	fwrite(buffer,sizeof(char),80,fp);

	// write geometry once
	if (outputType_ == GNUPLOT) text.open(geom_file_text.c_str(),ios_base::out);
	}
	else {
	fp=fopen(geom_file_name.c_str(),"ab"); // append

	// if (outputType_ == GNUPLOT) text.open(geom_file_text.c_str(),ios_base::app);
	}
	if (fp == NULL) {
	throw ATC_Error(0,"can not create Ensight geometry file");
	}

	// write preamble
	strcpy(buffer,"BEGIN TIME STEP");
	fwrite(buffer,sizeof(char),80,fp);
	strcpy(buffer,"Ensight geometry file");
	fwrite(buffer,sizeof(char),80,fp);
	strcpy(buffer,"description");
	fwrite(buffer,sizeof(char),80,fp);
	strcpy(buffer,"node id assign");
	fwrite(buffer,sizeof(char),80,fp);
	strcpy(buffer,"element id assign");
	fwrite(buffer,sizeof(char),80,fp);

	// per part
	strcpy(buffer,"part");
	fwrite(buffer,sizeof(char),80,fp);
	int part_number=1;
	fwrite(&part_number,sizeof(int),1,fp);
	strcpy(buffer,"description");
	fwrite(buffer,sizeof(char),80,fp);

	// write coordinates
	strcpy(buffer,"coordinates");
	fwrite(buffer,sizeof(char),80,fp);
	number_of_nodes_ = coordinates.nCols();
	fwrite(&number_of_nodes_,sizeof(int),1,fp);
	int number_of_spatial_dimensions = coordinates.nRows();
	if (number_of_spatial_dimensions != 3)
	throw ATC_Error(0,"Ensight writer needs a 3D geometry");

	for (int i = 0; i < number_of_spatial_dimensions; ++i)
	{
	for (int j = 0; j < number_of_nodes_; ++j)
	{
	float x = (float) coordinates(i,j);
	fwrite(&x,sizeof(float),1,fp);
	}
	}

	// write mesh connectivities or point "connectivities"
	if (connectivities)
	{
	dataType_ = MESH;
	strcpy(buffer,"hexa8");
	fwrite(buffer,sizeof(char),80,fp);
	int number_of_elements = connectivities->nCols();
	fwrite(&number_of_elements,sizeof(int),1,fp);
	int number_of_nodes_per_element = connectivities->nRows();
	for (int j = 0; j < number_of_elements; ++j)
	{
	for (int i = 0; i < number_of_nodes_per_element; ++i)
	{
	int inode = (*connectivities)(i,j) +1; // 1 based numbering
	fwrite(&inode,sizeof(int),1,fp);
	}
	}
	}
	else
	{
	strcpy(buffer,"point");
	fwrite(buffer,sizeof(char),80,fp);
	int number_of_elements = number_of_nodes_;
	fwrite(&number_of_elements,sizeof(int),1,fp);
	int number_of_nodes_per_element = 1;
	for (int j = 0; j < number_of_elements; ++j)
	{
	int inode = j +1; // 1 based numbering
	fwrite(&inode,sizeof(int),1,fp);
	}
	}

	// end per part
	strcpy(buffer,"END TIME STEP");
	fwrite(buffer,sizeof(char),80,fp);
	fclose(fp);

	// text output
	if (outputType_ == GNUPLOT )
	{
	for (int j = 0; j < number_of_nodes_; ++j)
	{
	for (int i = 0; i < number_of_spatial_dimensions; ++i)
	{
	text << setw(kFieldWidth) << std::scientific << std::setprecision(kFieldPrecison) << coordinates(i,j) << " ";
	}
	text << "\n";
	}
	text << "\n";
	}

	if (!initialized_) initialized_ = true;
	}
	//-----------------------------------------------------------------------------
	//*
	//-----------------------------------------------------------------------------
	void OutputManager::write_geometry(OUTPUT_LIST &part_coordinates)
	{
	if ( outputPrefix_ == "NULL") return;
	// geometry based on a reference configuration
	string geom_file_name = outputPrefix_ + ".geo";
	string geom_file_text = outputPrefix_ + ".XYZ";

	// open file
	FILE * fp =NULL;
	ofstream text;
	char buffer[80];
	if ( ! initialized_ )
	{
	fp=fopen(geom_file_name.c_str(),"wb"); // open
	strcpy(buffer,"C Binary");
	fwrite(buffer,sizeof(char),80,fp);

	if (outputType_ == GNUPLOT) text.open(geom_file_text.c_str(),ios_base::out);
	}
	else
	{
	fp=fopen(geom_file_name.c_str(),"ab"); // append
	if (outputType_ == GNUPLOT) text.open(geom_file_text.c_str(),ios_base::app);
	}
	if (fp == NULL) {
	throw ATC_Error(0,"can not create Ensight geometry file");
	}

	// write preamble
	strcpy(buffer,"BEGIN TIME STEP");
	fwrite(buffer,sizeof(char),80,fp);
	strcpy(buffer,"Ensight geometry file");
	fwrite(buffer,sizeof(char),80,fp);
	strcpy(buffer,"description");
	fwrite(buffer,sizeof(char),80,fp);
	strcpy(buffer,"node id assign");
	fwrite(buffer,sizeof(char),80,fp);
	strcpy(buffer,"element id assign");
	fwrite(buffer,sizeof(char),80,fp);

	// per part
	int part_number=1;
	OUTPUT_LIST::iterator iter;
	for (iter = part_coordinates.begin(); iter != part_coordinates.end(); iter++)
	{
	string part_name = iter->first;
	DenseMatrix<double> coordinates = *(iter->second);
	strcpy(buffer,"part");
	fwrite(buffer,sizeof(char),80,fp);
	fwrite(&part_number,sizeof(int),1,fp);
	strcpy(buffer,part_name.c_str());
	fwrite(buffer,sizeof(char),80,fp);

	// write coordinates
	strcpy(buffer,"coordinates");
	fwrite(buffer,sizeof(char),80,fp);
	number_of_nodes_ = coordinates.nCols();
	fwrite(&number_of_nodes_,sizeof(int),1,fp);
	int number_of_spatial_dimensions = coordinates.nRows();
	if (number_of_spatial_dimensions != 3)
	{
	throw ATC_Error(0,"Ensight writer needs a 3D geometry");
	}
	for (int i=0; i<number_of_spatial_dimensions; i++)
	{
	for (int j=0; j<number_of_nodes_; j++)
	{
	float x = (float)coordinates(i,j);
	fwrite(&x,sizeof(float),1,fp);
	}
	}

	// write mesh connectivities or point "connectivities"
	strcpy(buffer,"point");
	fwrite(buffer,sizeof(char),80,fp);
	int number_of_elements = number_of_nodes_;
	fwrite(&number_of_elements,sizeof(int),1,fp);
	int number_of_nodes_per_element = 1;
	for (int j = 0; j < number_of_elements; ++j)
	{
	int inode = j +1; // 1 based numbering
	fwrite(&inode,sizeof(int),1,fp);
	}
	++part_number;
	}

	strcpy(buffer,"END TIME STEP");
	fwrite(buffer,sizeof(char),80,fp);
	fclose(fp);

	if (!initialized_) initialized_ = true;
	}
	//-----------------------------------------------------------------------------
	/** pack "soln" into data */
	//-----------------------------------------------------------------------------
	void OutputManager::write_data(double time, FIELDS soln, OUTPUT_LIST data, const int *node_map)
	{
	// pack
	OUTPUT_LIST combined_data;
	if (soln)
	{
	FIELDS::iterator iter;
	for (iter = soln->begin(); iter != soln->end(); iter++)
	{
	FieldName field_index = iter->first;
	MATRIX* field_data = &(iter->second);
	string field_name = field_to_string(field_index);
	combined_data[field_name] = field_data;
	}
	}
	if (data)
	{
	OUTPUT_LIST::iterator iter;
	for (iter = data->begin(); iter != data->end(); iter++)
	{
	string field_name = iter->first;
	MATRIX* field_data = iter->second;
	combined_data[field_name] = field_data;
	}
	}
	write_data(time, &(combined_data), node_map);
	};
	//-----------------------------------------------------------------------------
	/** write (ensight gold format "C" binary) data */
	//-----------------------------------------------------------------------------
	void OutputManager::write_data(double time, OUTPUT_LIST data, const int node_map)
	{
	if (! initialized_) throw ATC_Error(0,"must write geometry before data");

	// write data
	OUTPUT_LIST::iterator iter;
	for (iter = data->begin(); iter != data->end(); iter++)
	{
	string field_name = iter->first;
	const MATRIX* field_data = iter->second;
	write_data(field_name, field_data, node_map);
	}

	// write dictionary
	write_dictionary(time,data);

	// write text dump
	if (outputType_ == GNUPLOT) {
	write_text(data);
	if (firstStep_ && node_map) {
	string map_file_text = outputPrefix_ + ".MAP";
	ofstream text;
	text.open(map_file_text.c_str(),ios_base::out);
	for (int i=0; i< number_of_nodes_ ; i++) {
	text << node_map[i] << "\n";
	}
	text.close();
	}
	}

	// global variables
	if (! globalData_.empty()) write_globals();

	if (firstStep_) firstStep_ = false;
	}
	//-----------------------------------------------------------------------------
	/** write (ensight gold format "C" binary) data */
	//-----------------------------------------------------------------------------
	void OutputManager::write_data(string field_name, const MATRIX field_data, const int node_map)
	{
	int ndof = field_data->nCols();
	int col_start = 0;
	int col_end = ndof;
	string filenames[kFileNameSize];
	int nfiles = 1;
	filenames[0] = outputPrefix_ + "." + field_name;
	int type = data_type(ndof);
	if (use_component_names(type)){
	nfiles = ndof;
	string* component_names = get_component_names(type);
	for (int ifile = 0; ifile < nfiles; ++ifile)
	{
	string comp_name = field_name + component_names[ifile];
	filenames[ifile] = outputPrefix_ + "." + comp_name;
	}
	}

	for (int ifile = 0; ifile < nfiles; ++ifile)
	{
	// for vector/tensor to components
	if ( nfiles > 1 )
	{
	col_start = ifile;
	col_end = ifile+1;
	}

	// open or append data file
	string data_file_name = filenames[ifile];
	FILE * fp;
	if ( outputTimes_.size() == 0 ) {
	fp=fopen(data_file_name.c_str(),"wb"); // open
	}
	else {
	fp=fopen(data_file_name.c_str(),"ab"); // append
	}
	if (fp == NULL) {
	throw ATC_Error(0,"can not create Ensight data file: "+data_file_name);
	}

	// write data
	char buffer[80];
	strcpy(buffer,"BEGIN TIME STEP");
	fwrite(buffer,sizeof(char),80,fp);
	strcpy(buffer,"field name"); // NOTE could be the field name
	fwrite(buffer,sizeof(char),80,fp);

	// per part
	strcpy(buffer,"part");
	fwrite(buffer,sizeof(char),80,fp);
	int part_number = 1;
	fwrite(&part_number,sizeof(int),1,fp);
	strcpy(buffer,"coordinates");
	fwrite(buffer,sizeof(char),80,fp);
	if (node_map)
	{
	for (int j = col_start; j < col_end; ++j)
	{
	for (int i = 0; i < number_of_nodes_; ++i)
	{
	int inode = node_map[i];
	float x = (float) (*field_data)(inode,j);
	fwrite(&x,sizeof(float),1,fp);
	}
	}
	}
	else
	{
	for (int j = col_start; j < col_end; ++j)
	{
	for (int i = 0; i < field_data->nRows(); ++i)
	{
	float x = (float) (*field_data)(i,j);
	fwrite(&x,sizeof(float),1,fp);
	}
	}
	}
	// end per part

	strcpy(buffer,"END TIME STEP");
	fwrite(buffer,sizeof(char),80,fp);
	fclose(fp);
	}
	}

	void OutputManager::write_text(OUTPUT_LIST *data)
	{
	string data_file_text = outputPrefix_ + ".DATA";
	ofstream text;
	if (firstStep_) text.open(data_file_text.c_str(),ios_base::out);
	else text.open(data_file_text.c_str(),ios_base::app);

	// write data label header
	if (firstStep_)
	{
	text.width(6); text << "# index:1" << " "; // give an ordinate for gnuplot
	text.width(10); text << " step:2" << " ";
	int k =3;
	if (data)
	{
	OUTPUT_LIST::iterator iter;
	for (iter = data->begin(); iter != data->end(); iter++)
	{
	int ncols = iter->second->nCols();
	string field_name = iter->first;
	if (ncols == 1) {
	string str; stringstream out; out << ":" << k; str = out.str();
	field_name.append(str);
	text.width(kFieldWidth); text << field_name << " ";
	k++;
	}
	else {
	for (int i = 1; i <= ncols; i++) {
	string str; stringstream out; out <<"_"<<i<<":"<<k; str = out.str();
	string name = field_name;
	name.append(str);
	text.width(kFieldWidth); text << name << " ";
	k++;
	}
	}
	}
	} else { throw ATC_Error(0," data missing from output");}
	text << "\n";
	}
	text << "# timestep " << outputTimes_.size() << " : "
	<< outputTimes_[outputTimes_.size()-1] << "\n";

	int nrows = 0;
	OUTPUT_LIST::iterator iter;
	iter = data->begin();
	const MATRIX* field_data = iter->second;
	nrows = field_data->nRows();

	for (int i = 0; i < nrows; ++i)
	{
	text.width(6); text << i << " "; // give an ordinate for gnuplot
	text.width(10); text << outputTimes_.size() << " ";
	OUTPUT_LIST::iterator iter;
	for (iter = data->begin(); iter != data->end(); iter++)
	{
	const MATRIX* field_data = iter->second;
	for (int j = 0; j < field_data->nCols(); ++j)
	{
	text.width(kFieldWidth);
	text << setw(kFieldWidth) << std::scientific << std::setprecision(kFieldPrecison) << (*field_data)(i,j) << " ";
	}
	}
	text <<"\n";
	}
	text <<"\n";
	}

	/** write (ensight gold : ASCII "C" format) dictionary */
	void OutputManager::write_dictionary(double time, OUTPUT_LIST *data)
	{
	// store the time step value
	outputTimes_.push_back(time);

	// file names
	string dict_file_name = outputPrefix_ + ".case";
	string geom_file_name = outputPrefix_ + ".geo";

	// open file
	FILE * fp;
	if ((fp=fopen(dict_file_name.c_str(),"w")) == NULL)
	{
	throw ATC_Error(0,"can not create Ensight case file");
	}

	// write file
	fprintf(fp,"FORMAT\n");
	fprintf(fp,"type: ensight gold\n");
	fprintf(fp,"GEOMETRY\n");
	if ( dataType_ == POINT) {
	fprintf(fp,"model: 1 1 %s change_coords_only\n", geom_file_name.c_str());
	} else {
	fprintf(fp,"model: %s\n", geom_file_name.c_str());
	}
	fprintf(fp,"VARIABLE\n");

	// data types
	if (!data) throw ATC_Error(0,"no data for output");
	OUTPUT_LIST::iterator iter;
	int ncols = 0;
	for (iter = data->begin(); iter != data->end(); iter++) {
	string field_name = iter->first;
	string field_file = outputPrefix_ + "." + field_name;
	const MATRIX* field_data = iter->second;
	int fieldCols = field_data->nCols();
	ncols += fieldCols;
	int type = data_type(fieldCols);
	if (use_component_names(type)){
	string* component_names = get_component_names(type);
	for (int j = 0; j < fieldCols; ++j)
	{
	string comp_name = field_name + component_names[j];
	string comp_file = outputPrefix_ + "." + comp_name;
	fprintf(fp,"scalar per node: 1 1 %s %s\n",
	comp_name.c_str(),comp_file.c_str());
	}
	}
	else if (type == VECTOR_OUTPUT) {
	fprintf(fp,"vector per node: 1 1 %s %s\n",
	field_name.c_str(),field_file.c_str());
	}
	else if (type == SYM_TENSOR_OUTPUT) {
	fprintf(fp,"tensor symm per node: 1 1 %s %s\n",
	field_name.c_str(),field_file.c_str());
	}
	else if (type == TENSOR_OUTPUT) {
	fprintf(fp,"tensor asymm per node: 1 1 %s %s\n",
	field_name.c_str(),field_file.c_str());
	}
	else {
	fprintf(fp,"scalar per node: 1 1 %s %s\n",
	field_name.c_str(),field_file.c_str());
	}
	}

	if (!firstStep_ && ncols != nDataCols_) {
	throw ATC_Error(0,"number of columns of data has changed: start new output");
	}
	nDataCols_ = ncols;

	int nsteps = outputTimes_.size();
	fprintf(fp,"TIME\n");
	fprintf(fp,"time set: 1\n");
	fprintf(fp,"number of steps: %10d\n",nsteps);
	if ( dataType_ == POINT) {
	fprintf(fp,"filename start number: 0\n");
	fprintf(fp,"filename increment: 1\n");
	}
	fprintf(fp,"time values:\n");
	for (int j = 0; j < nsteps; ++j) {
	double t = outputTimes_[j];
	fprintf(fp,"%12.5e",t);
	if ((j+1)%6 == 0) fprintf(fp,"\n");
	}
	fprintf(fp,"\n");
	fprintf(fp,"FILE\n");
	fprintf(fp,"file set: 1\n");
	fprintf(fp,"number of steps: %10d\n",nsteps);
	fclose(fp);
	};

	} // end ATC namespace

OutputManager.cppNo OneTemporaryActions

File Metadata

OutputManager.cppView Options

Event Timeline

OutputManager.cpp
No OneTemporary
Actions

OutputManager.cpp
View Options