text:
#include <fcntl.h>
#include <unistd.h>

#include <algorithm>
#include <cmath>
#include <cstdlib>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <sstream>
#include <vector>

#include <boost/program_options.hpp>
#include <boost/test/unit_test.hpp>

#include <xolotl/interface/Interface.h>
#include <xolotl/perf/dummy/DummyTimer.h>
#include <xolotl/perf/os/OSTimer.h>
#include <xolotl/test/MPITestUtils.h>
#include <xolotl/test/SystemTestCase.h>
#include <xolotl/test/config.h>
#include <xolotl/util/Filesystem.h>

using xolotl::test::SystemTestOptions;
BOOST_GLOBAL_FIXTURE(SystemTestOptions);

namespace xolotl
{
namespace bpo = ::boost::program_options;

namespace test
{
SystemTestOptions::SystemTestOptions()
{
	auto argc = boost::unit_test::framework::master_test_suite().argc;
	auto argv = boost::unit_test::framework::master_test_suite().argv;

	bpo::options_description desc("System Test Options");
	desc.add_options()("help,h", "produce help message")("verbose,v",
		"show all standard output")("approve,a", "approve running test cases")(
		"time-all,t", "report xolotl run time for each test case");

	bpo::variables_map opts;
	bpo::store(bpo::parse_command_line(argc, argv, desc), opts);
	bpo::notify(opts);

	if (opts.count("help")) {
		std::cout << desc << std::endl;
		std::exit(0);
	}

	if (opts.count("verbose")) {
		SystemTestCase::_noRedirect = true;
	}

	if (opts.count("approve")) {
		SystemTestCase::_approve = true;
	}

	if (opts.count("time-all")) {
		SystemTestCase::_timeAll = true;
	}
}

class StdOutRedirect
{
public:
	StdOutRedirect() = default;

	StdOutRedirect(const std::string& outputFileName)
	{
		if (getMPIRank() == 0) {
			fflush(stdout);
			int fd1 = open(
				outputFileName.c_str(), O_WRONLY | O_CREAT | O_TRUNC, 0644);
			if (fd1 < 0) {
				throw std::runtime_error(
					"Failed to open " + outputFileName + " for writing");
			}
			_fileDesc = dup(STDOUT_FILENO);
			if (_fileDesc < 0) {
				throw std::runtime_error("Failed to duplicate standard output");
			}
			if (dup2(fd1, STDOUT_FILENO) < 0) {
				throw std::runtime_error("Failed to duplicate " +
					outputFileName + " to standard output");
			}
			close(fd1);

			_redirected = true;
		}
	}

	StdOutRedirect(StdOutRedirect&& from) :
		_redirected(from._redirected),
		_fileDesc(from._fileDesc)
	{
		from._redirected = false;
	}

	StdOutRedirect&
	operator=(StdOutRedirect&& from)
	{
		_redirected = from._redirected;
		_fileDesc = from._fileDesc;
		from._redirected = false;
		return *this;
	}

	~StdOutRedirect() noexcept(false)
	{
		if (getMPIRank() == 0) {
			if (_redirected) {
				fflush(stdout);
				if (dup2(_fileDesc, STDOUT_FILENO) < 0) {
					throw std::runtime_error(
						"Failed to reinstate standard output");
				}
			}
		}
	}

private:
	bool _redirected{false};
	int _fileDesc;
};

class ScopedTimer
{
public:
	ScopedTimer(const std::string& name, bool enable) :
		_name(name),
		_enable(enable)
	{
		if (_enable) {
			MPI_Barrier(MPI_COMM_WORLD);
			if (getMPIRank() == 0) {
				_timer = std::make_unique<perf::os::OSTimer>();
				_timer->start();
			}
		}
	}

	~ScopedTimer()
	{
		if (_enable) {
			MPI_Barrier(MPI_COMM_WORLD);
			if (_timer) {
				_timer->stop();
				std::cout << _name << ": ";
				std::cout << _timer->getValue() << 's' << std::endl;
			}
		}
	}

private:
	bool _enable;
	std::unique_ptr<perf::ITimer> _timer;
	std::string _name;
};

double
diff2Norm(
	const std::vector<double>& data, const std::vector<double>& expectedData)
{
	double diffNorm = 0.0;
	double expectNorm = 0.0;
	auto sz = data.size();
	for (std::size_t i = 0; i < sz; ++i) {
		auto diff = data[i] - expectedData[i];
		auto expect = expectedData[i];
		if (std::isnan(expectedData[i])) {
			expect = 0.0;
			if (std::isnan(data[i])) {
				diff = 0.0;
			}
		}
		if (std::isinf(expectedData[i])) {
			expect = 0.0;
			if (std::isinf(data[i])) {
				diff = 0.0;
			}
		}
		diffNorm += diff * diff;
		expectNorm += expect * expect;
	}
	diffNorm = std::sqrt(diffNorm);
	expectNorm = std::sqrt(expectNorm);
	diffNorm /= expectNorm;
	return diffNorm;
}

double
diffInfNorm(
	const std::vector<double>& data, const std::vector<double>& expectedData)
{
	double diffNorm = 0.0;
	double expectNorm = 0.0;
	auto sz = data.size();
	for (std::size_t i = 0; i < sz; ++i) {
		auto diff = std::abs(data[i] - expectedData[i]);
		auto expect = std::abs(expectedData[i]);
		if (std::isnan(expectedData[i])) {
			expect = 0.0;
			if (std::isnan(data[i])) {
				diff = 0.0;
			}
		}
		if (std::isinf(expectedData[i])) {
			expect = 0.0;
			if (std::isinf(data[i])) {
				diff = 0.0;
			}
		}
		diffNorm = std::max(diffNorm, diff);
		expectNorm = std::max(expectNorm, expect);
	}
	diffNorm /= expectNorm;
	return diffNorm;
}

double
computeDiffNorm(
	const std::vector<double>& data, const std::vector<double>& expectedData)
{
	return diff2Norm(data, expectedData);
}

std::vector<double>
readOutputFile(const std::string& fileName)
{
	std::ifstream ifs(fileName);
	if (!ifs) {
		throw std::runtime_error("Unable to open file: " + fileName);
	}

	std::string line, tmpStr;
	std::vector<double> ret;
	while (getline(ifs, line)) {
		std::stringstream(line) >> tmpStr;
		if (tmpStr[0] == '#') {
			continue;
		}

		std::stringstream ss(line);
		while (ss >> tmpStr) {
			auto tmp = atof(tmpStr.c_str());
			ret.push_back(tmp);
		}
	}

	return ret;
}

const std::string SystemTestCase::_dataDir = TO_STRING(XOLOTL_TEST_DATA_DIR);
const std::string SystemTestCase::_binDir = TO_STRING(XOLOTL_BUILD_DIR);
const std::string SystemTestCase::_defaultOutputFileName = "retentionOut.txt";

bool SystemTestCase::_noRedirect = false;
bool SystemTestCase::_approve = false;
bool SystemTestCase::_timeAll = false;

SystemTestCase::SystemTestCase(
	const std::string& caseName, const std::string& outputFileName) :
	_caseName(caseName),
	_outputFileName(outputFileName)
{
	if (_timeAll) {
		_enableTimer = true;
	}
}

SystemTestCase::SystemTestCase(const std::string& caseName) :
	SystemTestCase(caseName, _defaultOutputFileName)
{
}

bool
SystemTestCase::runXolotl() const
{
	// Redirect console output
	StdOutRedirect redir;
	if (!_noRedirect) {
		redir = StdOutRedirect{
			_binDir + "/test/system/stdout_" + _caseName + ".txt"};
	}

	// Construct command-line
	auto exec = _binDir + "/xolotl/xolotl";
	auto paramsFileName = _dataDir + "/params_" + _caseName + ".txt";
	int argc = 2;
	const char* argv[] = {exec.data(), paramsFileName.data()};
	try {
		xolotl::interface::XolotlInterface{argc, argv}.solveXolotl();
	}
	catch (const std::exception& e) {
		std::cout << e.what() << std::endl;
		return false;
	}
	catch (...) {
		return false;
	}

	return true;
}

void
SystemTestCase::checkOutput(const std::string& outputFileName,
	const std::string& expectedOutputFileName) const
{
	auto expectedData = readOutputFile(expectedOutputFileName);
	auto data = readOutputFile(outputFileName);
	BOOST_REQUIRE(expectedData.size() == data.size());
	auto diffNorm = computeDiffNorm(data, expectedData);
	BOOST_REQUIRE_SMALL(diffNorm, _tolerance);
}

void
SystemTestCase::run() const
{
	{
		ScopedTimer timer{_caseName, _enableTimer};
		BOOST_REQUIRE(runXolotl());
	}

	if (getMPIRank() == 0) {
		if (_approve) {
			xolotl::fs::copy_file("./" + _outputFileName,
				_dataDir + "/output/" + _caseName + ".txt",
				xolotl::fs::copy_option::overwrite_if_exists);
		}
		else {
			checkOutput("./" + _outputFileName,
				_dataDir + "/output/" + _caseName + ".txt");
		}
	}
}

void
SystemTestCase::copyFile(const std::string& fileName)
{
	xolotl::fs::copy_file(_dataDir + "/" + fileName,
		xolotl::fs::absolute(fileName),
		xolotl::fs::copy_option::overwrite_if_exists);
}
} // namespace test
} // namespace xolotl

max_stars_repo_name: ORNL-Fusion/xolotl

max_issues_repo_name: ORNL-Fusion/xolotl

max_forks_repo_name: ORNL-Fusion/xolotl


text:
#include "car.h"
#include "genetic.h"
#include "layer.h"

#include <iostream>
#include <Eigen/Geometry>

Car::Car(): m_rotationToOriginal(0.0), m_mapSet(false), m_droveDistance(0.0),
    m_formerDistance(0.0), m_accumulatedRotation(0.0), m_carSize{4}
{
    setSpeed( 0.0 );
    setPosition( Eigen::Vector2d(0.0, 0.0));
    setAcceleration(0.0);
    setDirection(Eigen::Vector2d(1.0, 0.0));
    setRotationSpeed(0.0);

    setMeasureAngles( {-80, -50.0, -15.0, 0.0, 15.0, 50.0, 80} );

    std::vector<unsigned int> map = {8,4,2};
    m_network = NetworkPtr( new Network(map) );

    m_killer.start();
}

Car::~Car()
{

}

double Car::getSpeed() const
{
    return m_speed;
}

void Car::setSpeed(double speed)
{
    m_speed = speed;
}

const Eigen::Vector2d &Car::getPosition() const
{
    return m_position;
}

void Car::setPosition(const Eigen::Vector2d &position)
{
    m_position = position;
}

double Car::getAcceleration() const
{
    return m_acceleration;
}

void Car::setAcceleration(double acceleration)
{
    m_acceleration = acceleration;
}

const Eigen::Vector2d &Car::getDirection() const
{
    return m_direction;
}

void Car::setDirection(const Eigen::Vector2d &direction)
{
    m_direction = direction;
    m_direction.normalize();
    m_rotationToOriginal = computeAngleBetweenVectors(Eigen::Vector2d(1.0,0.0), m_direction);
}

void Car::update()
{
    double animTime = getTimeSinceLastUpdate();

    // rotate first
    double thisRotation =  m_rotationSpeedRad * animTime;
    Eigen::Rotation2D<double> r( thisRotation );
    m_direction = r.toRotationMatrix() * m_direction;
    Eigen::Vector2d(1.0,0.0);
    m_rotationToOriginal = computeAngleBetweenVectors(Eigen::Vector2d(1.0,0.0), m_direction);

    // adjust speed
    double newSpeed = std::max(m_speed + animTime*getAcceleration(), 0.0);
    newSpeed = std::min(newSpeed,600.0);
    Eigen::Vector2d effectiveSpeed = m_direction * (newSpeed + m_speed) * 0.5;

    // set new position
    Eigen::Vector2d newPosition = m_position + animTime * effectiveSpeed;
    newPosition = handleCollision(m_position, newPosition);

    // adjust drove distance and accumulated rotation
    m_droveDistance += (newPosition - getPosition()).norm();
    m_accumulatedRotation += std::abs(thisRotation);

    // important: measure distances before navigate and post move collision
    m_measuredDistances = measureDistances();
    navigate();

    // update
    m_speed = newSpeed;
    m_position = newPosition;
}

double Car::getFitness()
{
    return m_droveDistance;
}

double Car::getRotationSpeed() const
{
    return m_rotationSpeed;
}

void Car::setRotationSpeed(double rotationSpeed)
{
    m_rotationSpeed = rotationSpeed;
    m_rotationSpeedRad = -(M_PI / 180.0 * m_rotationSpeed);
}

double Car::getRotationRelativeToInitial() const
{
    return m_rotationToOriginal;
}

double Car::computeAngleBetweenVectors(const Eigen::Vector2d &a, const Eigen::Vector2d &b) const
{
    return std::acos( a.dot(b) / ( a.norm() * b.norm() ) ) / M_PI * 180.0;
}

std::shared_ptr<TrackMap> Car::getMap() const
{
    return m_map;
}

void Car::setMap(std::shared_ptr<TrackMap> map)
{
    m_map = map;
    m_mapSet = true;
}

Eigen::Vector2d Car::handleCollision(const Eigen::Vector2d& from, const Eigen::Vector2d& to)
{
    if( !m_mapSet )
        return to;

    if( !( m_map->isPositionValid(from) ) )
    {
        m_alive = false;
        return from;
    }

    Eigen::Vector2d dif = to - from;
    double l = dif.norm();

    if( l < 0.00000001 ) // no or very small move
        return to;

    Eigen::Vector2d du = dif * 1.0 / l;

    double tillEdge = distanceToEdge(from,du);

    if( tillEdge < l + m_carSize ) // collision
    {
        m_alive = false;
        return from;
    }

    return to;
}

void Car::handlePostMoveCollision()
{
    for(size_t k = 0; k < m_measuredDistances.rows(); k++)
    {
        if( m_measuredDistances(k,0) < m_carSize )
        {
            this->kill();
            return;
        }
    }
}

double Car::distanceToEdge(const Eigen::Vector2d &pos, const Eigen::Vector2d &direction) const
{
    if( !m_mapSet )
        return 0.0;

    Eigen::Vector2d d = direction.normalized();
    Eigen::Vector2d end = pos;

    bool goOn = true;

    int minDistToEdge = m_map->isPositionValid(end);

    while(minDistToEdge  > 0)
    {
        end = end + (d * minDistToEdge);
        minDistToEdge = m_map->isPositionValid(end);
    }

    return (end-pos).norm();
}

const std::vector<double> &Car::getMeasureAngles() const
{
    return m_measureAngles;
}

void Car::setMeasureAngles(const std::vector<double> &measureAngles)
{
    m_measureAngles = measureAngles;
}

Eigen::MatrixXd Car::measureDistances() const
{
    // distance, xP, yP
    Eigen::MatrixXd angs( m_measureAngles.size(), 3);
    for( size_t k = 0; k < m_measureAngles.size(); k++ )
    {
        double aRad = -(M_PI / 180.0 * m_measureAngles[k]);
        Eigen::Rotation2D<double> r( aRad );
        Eigen::Vector2d measureDir = r.toRotationMatrix() * m_direction;

        double dist = distanceToEdge(getPosition(),measureDir);

        Eigen::Vector2d edgePos = getPosition() + measureDir * dist;
        angs(k,0) = dist;
        angs(k,1) = edgePos(0);
        angs(k,2) = edgePos(1);
    }

    return angs;
}

Eigen::MatrixXd Car::getMeasuredDistances() const
{
    return m_measuredDistances;
}

void Car::considerSuicide()
{
    if(( getAge() > 1.0 && m_droveDistance < 3.0) || m_droveDistance < m_formerDistance * 1.01 )
    {
        m_alive = false;
    }
    else
    {
        m_formerDistance = m_droveDistance;
    }

    if( m_accumulatedRotation / getAge() > M_PI / 2.0) //90° per seconds average
    {
        m_alive = false;
    }
}

void Car::navigate()
{
    // decide what to do next
    Eigen::MatrixXd nnInput = m_measuredDistances.col(0);
    nnInput.conservativeResize(m_measuredDistances.rows()+1, 1); // additional input for speed
    nnInput(nnInput.rows()-1,0) = m_speed;

    // normalize input -> all values are positive -> scale them on a range -1 to +1
    double maxValInput = nnInput.maxCoeff();
    nnInput = (nnInput * 2.0/maxValInput).array() - 1.0;

    m_network->feedForward(nnInput);
    Eigen::MatrixXd nnOut = m_network->getOutputActivation();

    double maxRotationSpeed = 720.0;
    double maxAcceleration = 100.0;

    // scale output from 0 - 1 to -1 to +1
    double speedActivation = (nnOut(0,0) - 0.5) * 2;
    double rotationActivation = (nnOut(1,0) - 0.5) * 2;
    setAcceleration(maxAcceleration*speedActivation);
    setRotationSpeed(maxRotationSpeed*rotationActivation);

    if( m_killer.elapsed() > 1000 )
    {
        considerSuicide();
        m_killer.restart();
    }
}

ext: cpp

max_stars_repo_name: eidelen/EidNN

max_issues_repo_name: eidelen/EidNN

max_forks_repo_name: eidelen/EidNN


text:
#pragma once

#include <boost/mpl/identity.hpp>
#include <deal.II/base/graph_coloring.h>
#include <deal.II/dofs/dof_handler.h>

namespace boltzmann {
namespace local_ {

template <typename DH>
struct make_graph_coloring_traits
{
};

template <int DIM>
struct make_graph_coloring_traits<dealii::DoFHandler<DIM> >
{
  typedef typename dealii::DoFHandler<DIM>::active_cell_iterator active_iterator_t;
};
}  // end namespace _local

template <typename DH, typename INDEXER>
std::vector<std::vector<typename local_::make_graph_coloring_traits<DH>::active_iterator_t> >
make_graph_coloring(const DH& dof_handler, const INDEXER& indexer)
{
  // number of dofs...
  const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;

  auto begin = dof_handler.begin_active();
  auto end = dof_handler.end();

  typedef decltype(begin) iterator_t;

  /*
   *   The function make_graph_coloring creates a colored graph of the FEM Matrix
   *
   *   Hint: indexer.to_global is used to take into account periodic boundary
   *   conditions (if present).
   */
  const std::function<std::vector<dealii::types::global_dof_index>(const iterator_t&)>
      get_conflict_indices = [&](const iterator_t& cell) {
        std::vector<dealii::types::global_dof_index> dof_idxs(dofs_per_cell);
        // apply indexer, set velocity index j = 0;
        cell->get_dof_indices(dof_idxs);

        const unsigned int j = 0;
        for (unsigned int i = 0; i < dofs_per_cell; ++i) {
          dof_idxs[i] = indexer.to_global(dof_idxs[i], j);
        }
        return dof_idxs;
      };

  return dealii::GraphColoring::make_graph_coloring(begin, end, get_conflict_indices);
}

}  // end namespace boltzmann

ext: hpp


text:
//
// Copyright 2013 Christian Henning
//
// Distributed under the Boost Software License, Version 1.0
// See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt
//
#include <boost/gil.hpp>
#include <boost/gil/extension/io/jpeg.hpp>

#include <boost/core/lightweight_test.hpp>

#include "paths.hpp"
#include "scanline_read_test.hpp"

#include <cmath>
#include <sstream>
#include <string>

namespace gil = boost::gil;

#ifdef BOOST_GIL_IO_TEST_ALLOW_READING_IMAGES

void test_read_header() {
  using backend_t =
      gil::get_reader_backend<std::string const, gil::jpeg_tag>::type;
  backend_t backend = gil::read_image_info(jpeg_filename, gil::jpeg_tag());

  BOOST_TEST_EQ(backend._info._width, 1000u);
  BOOST_TEST_EQ(backend._info._height, 600u);

  BOOST_TEST_EQ(backend._info._num_components, 3);
  BOOST_TEST_EQ(backend._info._color_space, JCS_YCbCr);

  BOOST_TEST_EQ(backend._info._data_precision, 8);
}

void test_read_pixel_density() {
  using backend_t =
      gil::get_reader_backend<std::string const, gil::jpeg_tag>::type;
  backend_t backend =
      gil::read_image_info(jpeg_in + "EddDawson/36dpi.jpg", gil::jpeg_tag());

  gil::rgb8_image_t img;
  gil::read_image(jpeg_in + "EddDawson/36dpi.jpg", img, gil::jpeg_tag());

  gil::image_write_info<gil::jpeg_tag> write_settings;
  write_settings.set_pixel_dimensions(
      backend._info._width, backend._info._height,
      backend._info._pixel_width_mm, backend._info._pixel_height_mm);

  std::stringstream in_memory(std::ios_base::in | std::ios_base::out |
                              std::ios_base::binary);
  gil::write_view(in_memory, gil::view(img), write_settings);

  using backend2_t =
      gil::get_reader_backend<std::stringstream, gil::jpeg_tag>::type;
  backend2_t backend2 = gil::read_image_info(in_memory, gil::jpeg_tag());

  // Because of rounding the two results differ slightly.
  if (std::abs(backend._info._pixel_width_mm - backend2._info._pixel_width_mm) >
          10.0 ||
      std::abs(backend._info._pixel_height_mm -
               backend2._info._pixel_height_mm) > 10.0) {
    BOOST_TEST_EQ(0, 1);
  }
}

void test_read_reference_images() {
  using image_t = gil::rgb8_image_t;
  image_t img;
  read_image(jpeg_filename, img, gil::jpeg_tag());

#ifdef BOOST_GIL_IO_TEST_ALLOW_WRITING_IMAGES
  gil::write_view(jpeg_out + "rgb8_test.jpg", gil::view(img), gil::jpeg_tag());
#endif // BOOST_GIL_IO_TEST_ALLOW_WRITING_IMAGES
}

void test_dct_method_read() {
  using image_t = gil::rgb8_image_t;
  image_t img;

  gil::image_read_settings<gil::jpeg_tag> settings;
  settings._dct_method = gil::jpeg_dct_method::fast;

  gil::read_image(jpeg_filename, img, settings);

#ifdef BOOST_GIL_IO_TEST_ALLOW_WRITING_IMAGES
  gil::write_view(jpeg_out + "fast_dct_read_test.jpg", gil::view(img),
                  gil::jpeg_tag());
#endif // BOOST_GIL_IO_TEST_ALLOW_WRITING_IMAGES
}

void test_read_reference_images_image_iterator() {
  test_scanline_reader<gil::rgb8_image_t, gil::jpeg_tag>(jpeg_filename.c_str());
}

int main() {
  test_read_header();
  test_read_pixel_density();
  test_read_reference_images();
  test_dct_method_read();
  test_read_reference_images_image_iterator();

  return boost::report_errors();
}

#else
int main() {}
#endif // BOOST_GIL_IO_TEST_ALLOW_READING_IMAGES

max_stars_repo_name: sdebionne/gil-reformated

max_issues_repo_name: sdebionne/gil-reformated

max_forks_repo_name: sdebionne/gil-reformated


text:
#include "constants.h"
#include "inset_state.h"
#include "interpolate_bilinearly.h"
#include <boost/multi_array.hpp>
#include <omp.h>

// Function to calculate the velocity at the grid points (x, y) with x =
// 0.5, 1.5, ..., lx-0.5 and y = 0.5, 1.5, ..., ly-0.5 at time t
void calculate_velocity(
    double t,
    FTReal2d &grid_fluxx_init,
    FTReal2d &grid_fluxy_init,
    FTReal2d &rho_ft,
    FTReal2d &rho_init,
    boost::multi_array<double, 2> *grid_vx,
    boost::multi_array<double, 2> *grid_vy,
    const unsigned int lx,
    const unsigned int ly
) {
  double rho;

#pragma omp parallel for private(rho)
  for (unsigned int i = 0; i < lx; ++i) {
    for (unsigned int j = 0; j < ly; ++j) {
      rho = rho_ft(0, 0) + (1.0 - t) * (rho_init(i, j) - rho_ft(0,0));
      (*grid_vx)[i][j] = -grid_fluxx_init(i, j) / rho;
      (*grid_vy)[i][j] = -grid_fluxy_init(i, j) / rho;
    }
  }
  return;
}

bool all_points_are_in_domain(
    double delta_t,
    boost::multi_array<XYPoint, 2> *proj,
    boost::multi_array<XYPoint, 2> *v_intp,
    const unsigned int lx,
    const unsigned int ly
) {
  // Return false if and only if there exists a point that would be outside
  // [0, lx] x [0, ly]
  for (unsigned int i = 0; i < lx; ++i) {
    for (unsigned int j = 0; j < ly; ++j) {
      double x = (*proj)[i][j].x + 0.5 * delta_t * (*v_intp)[i][j].x;
      double y = (*proj)[i][j].y + 0.5 * delta_t * (*v_intp)[i][j].y;
      if (x < 0.0 || x > lx || y < 0.0 || y > ly) {
        return false;
      }
    }
  }
  return true;
}

// Function to integrate the equations of motion with the fast flow-based
// method
void InsetState::flatten_density()
{
  std::cerr << "In flatten_density()" << std::endl;

  // Constants for the numerical integrator
  const double inc_after_acc = 1.1;
  const double dec_after_not_acc = 0.75;
  const double abs_tol = (std::min(lx_, ly_) * 1e-6);

  // Resize proj_ multi-array if running for the first time
  if (proj_.shape()[0] != lx_ || proj_.shape()[1] != ly_) {
    proj_.resize(boost::extents[lx_][ly_]);
  }
  for (unsigned int i = 0; i < lx_; ++i) {
    for (unsigned int j = 0; j < ly_; ++j) {
      proj_[i][j].x = i + 0.5;
      proj_[i][j].y = j + 0.5;
    }
  }

  // Allocate memory for the velocity grid
  boost::multi_array<double, 2> grid_vx(boost::extents[lx_][ly_]);
  boost::multi_array<double, 2> grid_vy(boost::extents[lx_][ly_]);

  // Prepare Fourier transforms for the flux
  FTReal2d grid_fluxx_init;
  FTReal2d grid_fluxy_init;
  grid_fluxx_init.allocate(lx_, ly_);
  grid_fluxy_init.allocate(lx_, ly_);
  grid_fluxx_init.make_fftw_plan(FFTW_RODFT01, FFTW_REDFT01);
  grid_fluxy_init.make_fftw_plan(FFTW_REDFT01, FFTW_RODFT01);

  // eul[i][j] will be the new position of proj_[i][j] proposed by a simple
  // Euler step: move a full time interval delta_t with the velocity at time t
  // and position (proj_[i][j].x, proj_[i][j].y)
  boost::multi_array<XYPoint, 2> eul(boost::extents[lx_][ly_]);

  // mid[i][j] will be the new displacement proposed by the midpoint
  // method (see comment below for the formula)
  boost::multi_array<XYPoint, 2> mid(boost::extents[lx_][ly_]);

  // (vx_intp, vy_intp) will be the velocity at position (proj_.x, proj_.y) at
  // time t
  boost::multi_array<XYPoint, 2> v_intp(boost::extents[lx_][ly_]);

  // (vx_intp_half, vy_intp_half) will be the velocity at the midpoint
  // (proj_.x + 0.5*delta_t*vx_intp, proj_.y + 0.5*delta_t*vy_intp) at time
  // t + 0.5*delta_t
  boost::multi_array<XYPoint, 2> v_intp_half(boost::extents[lx_][ly_]);

  // Initialize the Fourier transforms of gridvx[] and gridvy[] at
  // every point on the lx_-times-ly_ grid at t = 0. We must typecast lx_ and ly_
  // as double-precision numbers. Otherwise the ratios in the denominator
  // will evaluate as zero.
  double dlx = lx_;
  double dly = ly_;

  // We temporarily insert the Fourier coefficients for the x-components and
  // y-components of the flux vector into grid_fluxx_init and grid_fluxy_init.
  // The reason for `+1` in `di+1` stems from the RODFT10 formula at:
  // https://www.fftw.org/fftw3_doc/1d-Real_002dodd-DFTs-_0028DSTs_0029.html
  for (unsigned int i = 0; i < lx_-1; ++i) {
    double di = i;
    for (unsigned int j = 0; j < ly_; ++j) {
      double denom = pi * ((di+1)/dlx + (j/(di+1)) * (j/dly) * (dlx/dly));
      grid_fluxx_init(i, j) =
        -rho_ft_(i+1, j) / denom;
    }
  }
  for (unsigned int j = 0; j < ly_; ++j) {
    grid_fluxx_init(lx_-1, j) = 0.0;
  }
  for (unsigned int i=0; i<lx_; ++i) {
    double di = i;
    for (unsigned int j = 0; j < ly_-1; ++j) {
      double denom = pi * ((di/(j+1)) * (di/dlx) * (dly/dlx) + (j+1)/dly);
      grid_fluxy_init(i, j) =
        -rho_ft_(i, j+1) / denom;
    }
  }
  for (unsigned int i=0; i<lx_; ++i) {
    grid_fluxy_init(i, ly_-1) = 0.0;
  }

  // Compute the flux vector and store the result in grid_fluxx_init and
  // grid_fluxy_init
  grid_fluxx_init.execute_fftw_plan();
  grid_fluxy_init.execute_fftw_plan();
  double t = 0.0;
  double delta_t = 1e-2;  // Initial time step.
  unsigned int iter = 0;

  // Integrate
  while (t < 1.0) {
    calculate_velocity(
      t,
      grid_fluxx_init,
      grid_fluxy_init,
      rho_ft_,
      rho_init_,
      &grid_vx,
      &grid_vy,
      lx_,
      ly_
    );
#pragma omp parallel for
    for (unsigned int i = 0; i < lx_; ++i) {
      for (unsigned int j = 0; j < ly_; ++j) {

        // We know, either because of the initialization or because of the
        // check at the end of the last iteration, that (proj_.x, proj_.y)
        // is inside the rectangle [0, lx_] x [0, ly_]. This fact guarantees
        // that interpolate_bilinearly() is given a point that cannot cause it
        // to fail.
        v_intp[i][j].x = interpolate_bilinearly(
          proj_[i][j].x,
          proj_[i][j].y,
          &grid_vx,
          'x',
          lx_,
          ly_
        );
        v_intp[i][j].y = interpolate_bilinearly(
          proj_[i][j].x,
          proj_[i][j].y,
          &grid_vy,
          'y',
          lx_,
          ly_
        );
      }
    }
    bool accept = false;
    while (!accept) {

      // Simple Euler step.

#pragma omp parallel for
      for (unsigned int i = 0; i < lx_; ++i) {
        for (unsigned int j = 0; j < ly_; ++j) {
          eul[i][j].x = proj_[i][j].x + v_intp[i][j].x * delta_t;
          eul[i][j].y = proj_[i][j].y + v_intp[i][j].y * delta_t;
        }
      }

      // Use "explicit midpoint method"
      // x <- x + delta_t * v_x(x + 0.5*delta_t*v_x(x,y,t),
      //                        y + 0.5*delta_t*v_y(x,y,t),
      //                        t + 0.5*delta_t)
      // and similarly for y.
      calculate_velocity(
        t + 0.5*delta_t,
        grid_fluxx_init,
        grid_fluxy_init,
        rho_ft_,
        rho_init_,
        &grid_vx,
        &grid_vy,
        lx_,
        ly_
      );

      // Make sure we do not pass a point outside [0, lx_] x [0, ly_] to
      // interpolate_bilinearly(). Otherwise decrease the time step below and
      // try again.
      accept = all_points_are_in_domain(delta_t, &proj_, &v_intp, lx_, ly_);
      if (accept) {

        // Okay, we can run interpolate_bilinearly()

#pragma omp parallel for
        for (unsigned int i = 0; i < lx_; ++i) {
          for (unsigned int j = 0; j < ly_; ++j) {
            v_intp_half[i][j].x = interpolate_bilinearly(
              proj_[i][j].x + 0.5*delta_t*v_intp[i][j].x,
              proj_[i][j].y + 0.5*delta_t*v_intp[i][j].y,
              &grid_vx,
              'x',
              lx_,
              ly_
            );
            v_intp_half[i][j].y = interpolate_bilinearly(
              proj_[i][j].x + 0.5*delta_t*v_intp[i][j].x,
              proj_[i][j].y + 0.5*delta_t*v_intp[i][j].y,
              &grid_vy,
              'y',
              lx_,
              ly_
            );
            mid[i][j].x = proj_[i][j].x + v_intp_half[i][j].x * delta_t;
            mid[i][j].y = proj_[i][j].y + v_intp_half[i][j].y * delta_t;

            // Do not accept the integration step if the maximum squared
            // difference between the Euler and midpoint proposals exceeds
            // abs_tol. Neither should we accept the integration step if one
            // of the positions wandered out of the domain. If one of these
            // problems occurred, decrease the time step.
            const double sq_dist =
              (mid[i][j].x-eul[i][j].x) * (mid[i][j].x-eul[i][j].x)
              + (mid[i][j].y-eul[i][j].y) * (mid[i][j].y-eul[i][j].y);
            if (sq_dist > abs_tol ||
                mid[i][j].x < 0.0 || mid[i][j].x > lx_ ||
                mid[i][j].y < 0.0 || mid[i][j].y > ly_) {
              accept = false;
            }
          }
        }
      }
      if (!accept) {
        delta_t *= dec_after_not_acc;
      }
    }

    // Control ouput
    if (iter % 10 == 0) {
      std::cerr << "iter = "
                << iter
                << ", t = "
                << t
                << ", delta_t = "
                << delta_t
                << "\n";
    }

    // When we get here, the integration step was accepted
    t += delta_t;
    ++iter;
    proj_ = mid;
    delta_t *= inc_after_acc;  // Try a larger step next time
  }
  grid_fluxx_init.destroy_fftw_plan();
  grid_fluxy_init.destroy_fftw_plan();
  grid_fluxx_init.free();
  grid_fluxy_init.free();
  return;
}

max_stars_repo_name: mgastner/cartogram-cpp

max_issues_repo_name: mgastner/cartogram-cpp

max_forks_repo_name: mgastner/cartogram-cpp


text:
#ifndef UTIL_THREAD_POOL_H
#define UTIL_THREAD_POOL_H

#include "util/pcqueue.hh"

#include <boost/ptr_container/ptr_vector.hpp>
#include <boost/optional.hpp>
#include <boost/thread.hpp>

#include <iostream>
#include <cstdlib>

namespace util {

template <class HandlerT> class Worker : boost::noncopyable {
  public:
    typedef HandlerT Handler;
    typedef typename Handler::Request Request;

    template <class Construct> Worker(PCQueue<Request> &in, Construct &construct, const Request &poison)
      : in_(in), handler_(construct), poison_(poison), thread_(boost::ref(*this)) {}

    // Only call from thread.
    void operator()() {
      Request request;
      while (1) {
        in_.Consume(request);
        if (request == poison_) return;
        try {
          (*handler_)(request);
        }
        catch(const std::exception &e) {
          std::cerr << "Handler threw " << e.what() << std::endl;
          abort();
        }
        catch(...) {
          std::cerr << "Handler threw an exception, dropping request" << std::endl;
          abort();
        }
      }
    }

    void Join() {
      thread_.join();
    }

  private:
    PCQueue<Request> &in_;

    boost::optional<Handler> handler_;

    const Request poison_;

    boost::thread thread_;
};

template <class HandlerT> class ThreadPool : boost::noncopyable {
  public:
    typedef HandlerT Handler;
    typedef typename Handler::Request Request;

    template <class Construct> ThreadPool(size_t queue_length, size_t workers, Construct handler_construct, Request poison) : in_(queue_length), poison_(poison) {
      for (size_t i = 0; i < workers; ++i) {
        workers_.push_back(new Worker<Handler>(in_, handler_construct, poison));
      }
    }

    ~ThreadPool() {
      for (size_t i = 0; i < workers_.size(); ++i) {
        Produce(poison_);
      }
      for (typename boost::ptr_vector<Worker<Handler> >::iterator i = workers_.begin(); i != workers_.end(); ++i) {
        i->Join();
      }
    }

    void Produce(const Request &request) {
      in_.Produce(request);
    }

    // For adding to the queue.
    PCQueue<Request> &In() { return in_; }

  private:
    PCQueue<Request> in_;

    boost::ptr_vector<Worker<Handler> > workers_;

    Request poison_;
};

} // namespace util

#endif // UTIL_THREAD_POOL_H

ext: hh

max_stars_repo_name: ugermann/MMT

max_issues_repo_name: Pangeamt/nectm

max_forks_repo_name: Pangeamt/nectm


text:
//==============================================================================
//         Copyright 2003 - 2011 LASMEA UMR 6602 CNRS/Univ. Clermont II
//         Copyright 2009 - 2011 LRI    UMR 8623 CNRS/Univ Paris Sud XI
//
//          Distributed under the Boost Software License, Version 1.0.
//                 See accompanying file LICENSE.txt or copy at
//                     http://www.boost.org/LICENSE_1_0.txt
//==============================================================================
#ifndef BOOST_SIMD_BITWISE_FUNCTIONS_SIMD_SSE_AVX_SHRI_HPP_INCLUDED
#define BOOST_SIMD_BITWISE_FUNCTIONS_SIMD_SSE_AVX_SHRI_HPP_INCLUDED
#ifdef BOOST_SIMD_HAS_AVX_SUPPORT
#include <boost/simd/bitwise/functions/shri.hpp>
#include <boost/dispatch/meta/as_integer.hpp>

namespace boost { namespace simd { namespace ext
{
  BOOST_SIMD_FUNCTOR_IMPLEMENTATION( boost::simd::tag::shri_, boost::simd::tag::avx_
                            , (A0)(A1)
                            , ((simd_<arithmetic_<A0>,boost::simd::tag::avx_>))
                              (scalar_< integer_<A1> >)
                            )
  {
    typedef A0 result_type;
    BOOST_SIMD_FUNCTOR_CALL(2)
    {
      typedef typename meta::scalar_of<A0>::type                       sctype;
      typedef typename dispatch::meta::as_integer<sctype>::type        sitype;
      typedef simd::native<sitype, boost::simd::tag::sse_ >           isvtype;
      typedef typename dispatch::meta::as_integer<A0>::type             itype;

      isvtype a00 = _mm256_extractf128_si256(simd::bitwise_cast<itype>(a0), 0);
      isvtype a01 = _mm256_extractf128_si256(simd::bitwise_cast<itype>(a0), 1);
      itype that = _mm256_insertf128_si256(that,boost::simd::shri( a00, a1), 0);
      return  simd::bitwise_cast<A0>(_mm256_insertf128_si256(that, boost::simd::shri(a01, a1), 1));
     }
  };
} } }

#endif
#endif


text:
//------------------------------------------------------------------------------
/*
    This file is part of rippled: https://github.com/ripple/rippled
    Copyright (c) 2012-2014 Ripple Labs Inc.

    Permission to use, copy, modify, and/or distribute this software for any
    purpose  with  or without fee is hereby granted, provided that the above
    copyright notice and this permission notice appear in all copies.

    THE  SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
    WITH  REGARD  TO  THIS  SOFTWARE  INCLUDING  ALL  IMPLIED  WARRANTIES  OF
    MERCHANTABILITY  AND  FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
    ANY  SPECIAL ,  DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
    WHATSOEVER  RESULTING  FROM  LOSS  OF USE, DATA OR PROFITS, WHETHER IN AN
    ACTION  OF  CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
    OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
*/
//==============================================================================

#include <ripple/app/main/Application.h>
#include <ripple/app/misc/Transaction.h>
#include <ripple/core/DatabaseCon.h>
#include <ripple/core/SociDB.h>
#include <ripple/net/RPCErr.h>
#include <ripple/protocol/ErrorCodes.h>
#include <ripple/protocol/jss.h>
#include <ripple/resource/Fees.h>
#include <ripple/rpc/Context.h>
#include <ripple/rpc/Role.h>
#include <boost/format.hpp>

namespace ripple {

// {
//   start: <index>
// }
Json::Value doTxHistory (RPC::JsonContext& context)
{
    context.loadType = Resource::feeMediumBurdenRPC;

    if (!context.params.isMember (jss::start))
        return rpcError (rpcINVALID_PARAMS);

    unsigned int startIndex = context.params[jss::start].asUInt ();

    if ((startIndex > 10000) &&  (! isUnlimited (context.role)))
        return rpcError (rpcNO_PERMISSION);

    Json::Value obj;
    Json::Value txs;

    obj[jss::index] = startIndex;

    std::string sql =
        boost::str (boost::format (
            "SELECT LedgerSeq, Status, RawTxn "
            "FROM Transactions ORDER BY LedgerSeq desc LIMIT %u,20;")
                    % startIndex);

    {
        auto db = context.app.getTxnDB ().checkoutDb ();

        boost::optional<std::uint64_t> ledgerSeq;
        boost::optional<std::string> status;
        soci::blob sociRawTxnBlob (*db);
        soci::indicator rti;
        Blob rawTxn;

        soci::statement st = (db->prepare << sql,
                              soci::into (ledgerSeq),
                              soci::into (status),
                              soci::into (sociRawTxnBlob, rti));

        st.execute ();
        while (st.fetch ())
        {
            if (soci::i_ok == rti)
                convert(sociRawTxnBlob, rawTxn);
            else
                rawTxn.clear ();

            if (auto trans = Transaction::transactionFromSQL (
                    ledgerSeq, status, rawTxn, context.app))
                txs.append (trans->getJson (JsonOptions::none));
        }
    }

    obj[jss::txs] = txs;

    return obj;
}

} // ripple

max_stars_repo_name: tlongwell-ripple/rippled

max_issues_repo_name: tlongwell-ripple/rippled

max_forks_repo_name: tlongwell-ripple/rippled


text:
#pragma once
#include <boost/python.hpp>
#include <boost/circular_buffer.hpp>

namespace gm
{
    namespace Python
    {
        class Output
        {
        typedef boost::circular_buffer<std::string> ContainerType;
        private:
            static ContainerType m_outputs;
            unsigned int m_bufferSize = 1000;
        public:
            void setBufferSize(unsigned int buffer)
            {
                this->m_bufferSize = buffer;
            }

            void write( std::string const& str )
            {
                if (m_outputs.capacity()<m_bufferSize)
                    m_outputs.resize(m_bufferSize);
                m_outputs.push_back(str);
            }

            static std::string GetOutput()
            {
                std::string ret;
                std::stringstream ss;
                for(auto it=m_outputs.begin(); it!=m_outputs.end(); it++)
                {
                    ss << *it;
                }
                m_outputs.clear();
                ret =  ss.str();
                return ret;
            }
        };

        Output::ContainerType Output::m_outputs;
    }
}

max_stars_repo_name: GraphMIC/GraphMIC

max_issues_repo_name: kevinlq/GraphMIC

max_forks_repo_name: kevinlq/GraphMIC


text:
#include <iostream>
#include <boost/bind.hpp>
#include <boost/smart_ptr/scoped_ptr.hpp>
#include <boost/algorithm/string/split.hpp>
#include <boost/algorithm/string/trim.hpp>
#include <boost/algorithm/string/case_conv.hpp>
#include <boost/algorithm/string/predicate.hpp>
#include <boost/type_traits/conditional.hpp>
#include "cancellable_streambuf.h"
#include "data_receiver.h"
#include "sample.h"
#include "socket_utils.h"

// a convention that applies when including portable_oarchive.h in multiple .cpp files.
// otherwise, the templates are instantiated in this file and sample.cpp which leads
// to errors like "multiple definition of `typeinfo name"
#define NO_EXPLICIT_TEMPLATE_INSTANTIATION
#include "portable_archive/portable_iarchive.hpp"

// === implementation of the data_receiver class ===

using namespace lslboost::algorithm;

namespace lsl {

/**
* Construct a new data receiver from an info connection.
* @param conn An inlet connection object.
* @param max_buflen Optionally the maximum amount of data to buffer (in seconds if there is a nominal sampling rate, otherwise x 100 in samples).
*					Recording applications want to use a fairly large buffer size here, while real-time applications want to only buffer as much as they need to perform their next calculation.
* @param max_chunklen Optionally the maximum size, in samples, at which chunks are transmitted (the default corresponds to the chunk sizes used by the sender).
*					  Recording applications can use a generous size here (leaving it to the network how to pack things), while real-time applications may want a finer (perhaps 1-sample) granularity.
*/
data_receiver::data_receiver(inlet_connection &conn, int max_buflen, int max_chunklen): conn_(conn), check_thread_start_(true), closing_stream_(false), connected_(false), sample_queue_(max_buflen),
	sample_factory_(new factory(conn.type_info().channel_format(),conn.type_info().channel_count(),conn.type_info().nominal_srate()?conn.type_info().nominal_srate()*api_config::get_instance()->inlet_buffer_reserve_ms()/1000:api_config::get_instance()->inlet_buffer_reserve_samples())), max_buflen_(max_buflen), max_chunklen_(max_chunklen)
{
	if (max_buflen < 0)
		throw std::invalid_argument("The max_buflen argument must not be smaller than 0.");
	if (max_chunklen < 0)
		throw std::invalid_argument("The max_chunklen argument must not be smaller than 0.");
	conn_.register_onlost(this,&connected_upd_);
}

/// Destructor. Stops the background activities.
data_receiver::~data_receiver() {
	try {
		conn_.unregister_onlost(this);
		if (data_thread_.joinable())
			data_thread_.join();
	}
	catch(std::exception &e) {
		std::cerr << "Unexpected error during destruction of a data_receiver: " << e.what() << std::endl;
	}
	catch(...) {
		std::cerr << "Severe error during data receiver shutdown." << std::endl;
	}
}

// === external access ===

/**
* Open a new data stream.
* All samples pushed in at the other end from this moment onwards will be queued and
* eventually be delivered in response to pull_sample() or pull_chunk() calls.
* A pull call without preceding open_stream serves as an implicit open_stream.
* @param timeout Optional timeout of the operation (default: no timeout).
* @throws timeout_error (if the timeout expires), or lost_error (if the stream source has been lost).
*/
void data_receiver::open_stream(double timeout) {
	closing_stream_ = false;
	lslboost::unique_lock<lslboost::mutex> lock(connected_mut_);
	if (!connection_completed()) {
		// start thread if not yet running
		if (check_thread_start_ && !data_thread_.joinable()) {
			data_thread_ = lslboost::thread(&data_receiver::data_thread,this);
			check_thread_start_ = false;
		}
		// wait until the connection attempt completes (or we time out)
		if (timeout >= FOREVER)
			connected_upd_.wait(lock, lslboost::bind(&data_receiver::connection_completed,this));
		else
			if (!connected_upd_.wait_for(lock, lslboost::chrono::duration<double>(timeout), lslboost::bind(&data_receiver::connection_completed,this)))
				throw timeout_error("The open_stream() operation timed out.");
	}
	if (conn_.lost())
		throw lost_error("The stream read by this inlet has been lost. To recover, you need to re-resolve the source and re-create the inlet.");
}

/**
* Close the current data stream.
* All samples still buffered or in flight will be dropped and the source will halt its buffering of data for this inlet.
* If an application stops being interested in data from a source (temporarily or not), it should call drop_stream() to not
* pressure the source outlet to buffer unnecessarily large amounts of data (perhaps even running out of memory).
*/
void data_receiver::close_stream() {
	check_thread_start_ = true;
	closing_stream_ = true;
	cancel_all_registered();
}

template <class T>
double data_receiver::pull_sample_typed(T *buffer, int buffer_elements, double timeout) {
	if (conn_.lost())
		throw lost_error("The stream read by this outlet has been lost. To recover, you need to "
						 "re-resolve the source and re-create the inlet.");
	// start data thread implicitly if necessary
	if (check_thread_start_ && !data_thread_.joinable()) {
		data_thread_ = lslboost::thread(&data_receiver::data_thread, this);
		check_thread_start_ = false;
	}
	// get the sample with timeout
	if (sample_p s = sample_queue_.pop_sample(timeout)) {
		if (buffer_elements != conn_.type_info().channel_count())
			throw std::range_error("The number of buffer elements provided does not match the "
								   "number of channels in the sample.");
		s->retrieve_typed(buffer);
		return s->timestamp;
	} else {
		if (conn_.lost())
			throw lost_error("The stream read by this inlet has been lost. To recover, you need to "
							 "re-resolve the source and re-create the inlet.");
		return 0.0;
	}
}

typedef lslboost::conditional<sizeof(long) == 8, int64_t, int32_t>::type long_type;
template <>
double data_receiver::pull_sample_typed(long *buffer, int buffer_elements, double timeout) {
	return pull_sample_typed((long_type *)buffer, buffer_elements, timeout);
}

template double data_receiver::pull_sample_typed<char>(char *, int, double);
template double data_receiver::pull_sample_typed<int16_t>(int16_t *, int, double);
template double data_receiver::pull_sample_typed<int32_t>(int32_t *, int, double);
template double data_receiver::pull_sample_typed<int64_t>(int64_t *, int, double);
template double data_receiver::pull_sample_typed<float>(float *, int, double);
template double data_receiver::pull_sample_typed<double>(double *, int, double);
template double data_receiver::pull_sample_typed<std::string>(std::string *, int, double);

/**
* Pull a sample from the inlet and read it into a pointer to raw data.
* No type checking or conversions are done (not recommended!). Do not use for variable-size/string-formatted streams.
* @param buffer A pointer to hold the resulting raw sample data.
* @param buffer_bytes The number of bytes allocated in the buffer. Note: it is the responsibility of the user to allocate enough memory.
* @param block If true, the function will block until a sample is available; otherwise, it will return 0.0 and no data.
* @return The capture time of the sample on the remote machine, or 0.0 if no new sample was available.
*		   To remap this time stamp to the local clock, add the value returned by .time_correction() to it.
*		   This is only necessary if the clocks of the source and destination machine are not synchronized to high enough precision.
*/
double data_receiver::pull_sample_untyped(void *buffer, int buffer_bytes, double timeout) {
	if (conn_.lost())
		throw lost_error("The stream read by this inlet has been lost. To recover, you need to re-resolve the source and re-create the inlet.");
	// start data thread implicitly if necessary
	if (check_thread_start_ && !data_thread_.joinable()) {
		data_thread_ = lslboost::thread(&data_receiver::data_thread,this);
		check_thread_start_ = false;
	}
	// get the sample with timeout
	if (sample_p s = sample_queue_.pop_sample(timeout)) {
		if (buffer_bytes != conn_.type_info().sample_bytes())
			throw std::range_error("The size of the provided buffer does not match the number of bytes in the sample.");
		s->retrieve_untyped(buffer);
		return s->timestamp;
	} else {
		if (conn_.lost())
			throw lost_error("The stream read by this inlet has been lost. To recover, you need to re-resolve the source and re-create the inlet.");
		return 0.0;
	}
}

// === internal processing ===

/// The data reader thread.
void data_receiver::data_thread() {
	conn_.acquire_watchdog();
	// ensure that the sample factory persists for the lifetime of this thread
	factory_p factory(sample_factory_);
	try {
		while (!conn_.lost() && !conn_.shutdown() && !closing_stream_) {
			try {
				// --- connection setup ---

				// make a new stream buffer and a stream on top of it
				cancellable_streambuf buffer;
				buffer.register_at(&conn_);
				buffer.register_at(this);
				std::iostream server_stream(&buffer);
				lslboost::scoped_ptr<eos::portable_iarchive> inarch;
				// connect to endpoint
				buffer.connect(conn_.get_tcp_endpoint());
				if (buffer.puberror())
					throw buffer.puberror();

				// --- protocol negotiation ---

				int use_byte_order = 0;				// which byte order we shall use (0=portable byte order)
				int data_protocol_version = 100;	// which protocol version we shall use for data transmission (100=version 1.00)
				bool suppress_subnormals = false;	// whether we shall suppress subnormal numbers

				// propose to use the highest protocol version supported by both parties
				int proposed_protocol_version = std::min(api_config::get_instance()->use_protocol_version(),conn_.type_info().version());
				if (proposed_protocol_version >= 110) {
					// request line LSL:streamfeed/[ProtocolVersion] [UID]\r\n
					server_stream << "LSL:streamfeed/" << proposed_protocol_version << " " << conn_.current_uid() << "\r\n";
					// transmit request parameters
					server_stream << "Native-Byte-Order: " << BOOST_BYTE_ORDER << "\r\n";
					server_stream << "Endian-Performance: " << std::floor(measure_endian_performance()) << "\r\n";
					server_stream << "Has-IEEE754-Floats: " << (format_ieee754[cft_float32] && format_ieee754[cft_double64]) << "\r\n";
					server_stream << "Supports-Subnormals: " << format_subnormal[conn_.type_info().channel_format()] << "\r\n";
					server_stream << "Value-Size: " << conn_.type_info().channel_bytes() << "\r\n"; // 0 for strings
					server_stream << "Data-Protocol-Version: " << proposed_protocol_version << "\r\n";
					server_stream << "Max-Buffer-Length: " << max_buflen_ << "\r\n";
					server_stream << "Max-Chunk-Length: " << max_chunklen_ << "\r\n";
					server_stream << "Hostname: " << conn_.type_info().hostname() << "\r\n";
					server_stream << "Source-Id: " << conn_.type_info().source_id() << "\r\n";
					server_stream << "Session-Id: " << conn_.type_info().session_id() << "\r\n";
					server_stream << "\r\n" << std::flush;

					// check server response line (LSL/[Version] [StatusCode] [Message])
					char buf[16384] = {0};
					if (!server_stream.getline(buf,sizeof(buf)))
						throw lost_error("Connection lost.");
					std::vector<std::string> parts; split(parts,buf,is_any_of(" \t"));
					if (parts.size() < 3 || !starts_with(parts[0],"LSL/"))
						throw std::runtime_error("Received a malformed response.");
					if (from_string<int>(parts[0].substr(4))/100 > api_config::get_instance()->use_protocol_version()/100)
						throw std::runtime_error("The other party's protocol version is too new for this client; please upgrade your LSL library.");
					int status_code = from_string<int>(parts[1]);
					if (status_code == 404)
						throw lost_error("The given address does not serve the resolved stream (likely outdated).");
					if (status_code >= 400)
						throw std::runtime_error("The other party sent an error: " + std::string(buf));
					if (status_code >= 300)
						throw lost_error("The other party requested a redirect.");

					// receive response parameters
					while (server_stream.getline(buf,sizeof(buf)) && (buf[0] != '\r')) {
						std::string hdrline(buf);
						std::size_t colon = hdrline.find_first_of(':');
						if (colon != std::string::npos) {
							// extract key & value
							std::string type = to_lower_copy(trim_copy(hdrline.substr(0,colon))), rest = to_lower_copy(trim_copy(hdrline.substr(colon+1)));
							// strip off comments
							std::size_t semicolon = rest.find_first_of(';');
							if (semicolon != std::string::npos)
								rest = rest.substr(0,semicolon);
							// get the header information
							if (type == "byte-order") {
								use_byte_order = from_string<int>(rest);
								if (use_byte_order==2134 && BOOST_BYTE_ORDER!=2134 && format_sizes[conn_.type_info().channel_format()]>=8)
									throw std::runtime_error("The byte order conversion requested by the other party is not supported.");
							}
							if (type == "suppress-subnormals")
								suppress_subnormals = from_string<bool>(rest);
							if (type == "uid" && rest != conn_.current_uid())
								throw lost_error("The received UID does not match the current connection's UID.");
							if (type == "data-protocol-version") {
								data_protocol_version = from_string<int>(rest);
								if (data_protocol_version > api_config::get_instance()->use_protocol_version())
									throw std::runtime_error("The protocol version requested by the other party is not supported by this client.");
							}
						}
					}
					if (!server_stream)
						throw lost_error("Server connection lost.");
				} else {
					// version 1.00: send request line and feed parameters
					server_stream << "LSL:streamfeed\r\n";
					server_stream << max_buflen_ << " " << max_chunklen_ << "\r\n" << std::flush;
				}

				if (data_protocol_version == 100) {
					// portable binary archive (parse archive header)
					inarch.reset(new eos::portable_iarchive(server_stream));
					// receive stream_info message from server
					std::string infomsg; *inarch >> infomsg;
					stream_info_impl info; info.from_shortinfo_message(infomsg);
					// confirm that the UID matches, otherwise reconnect
					if (info.uid() != conn_.current_uid())
						throw lost_error("The received UID does not match the current connection's UID.");
				}

				// --- format validation ---
				{
					// receive and parse two subsequent test-pattern samples and check if they are formatted as expected
					lslboost::scoped_ptr<sample> temp[4];
					for (int k=0; k<4; temp[k++].reset(factory::new_sample_unmanaged(conn_.type_info().channel_format(),conn_.type_info().channel_count(),0.0,false)));
					temp[0]->assign_test_pattern(4);
					if (data_protocol_version >= 110)
						temp[1]->load_streambuf(buffer, data_protocol_version, use_byte_order, suppress_subnormals);
					else
						*inarch >> *temp[1];
					temp[2]->assign_test_pattern(2);
					if (data_protocol_version >= 110)
						temp[3]->load_streambuf(buffer, data_protocol_version, use_byte_order, suppress_subnormals);
					else
						*inarch >> *temp[3];
					if (!(*temp[0].get() == *temp[1].get()) || !(*temp[2].get() == *temp[3].get()))
						throw std::runtime_error("The received test-pattern samples do not match the specification. The protocol formats are likely incompatible.");
				}

                // signal to accessor functions on other threads that the protocol negotiation has been successful,
                // so we're now connected (and remain to be even if we later recover silently)
                {
                    lslboost::lock_guard<lslboost::mutex> lock(connected_mut_);
                    connected_ = true;
                }
                connected_upd_.notify_all();

                // --- transmission loop ---

                double last_timestamp = 0.0;
                double srate = conn_.current_srate();
                for (int k=0;!conn_.lost() && !conn_.shutdown() && !closing_stream_;k++) {
                    // allocate and fetch a new sample
                    sample_p samp(factory->new_sample(0.0,false));
					if (data_protocol_version >= 110) samp->load_streambuf(buffer,data_protocol_version,use_byte_order,suppress_subnormals); else *inarch >> *samp;
                    // deduce timestamp if necessary
                    if (samp->timestamp == DEDUCED_TIMESTAMP) {
                        samp->timestamp = last_timestamp;
                        if (srate != IRREGULAR_RATE)
                            samp->timestamp += 1.0/srate;
                    }
                    last_timestamp = samp->timestamp;
                    // push it into the sample queue
                    sample_queue_.push_sample(samp);
                    // periodically update the last receive time to keep the watchdog happy
					if (srate<=16 || (k & 0xF) == 0)
						conn_.update_receive_time(lsl_clock());
                }
			}
			catch(error_code &) {
				// connection-level error: closed, reset, refused, etc.
				conn_.try_recover_from_error();
			}
			catch(lost_error &) {
				// another type of connection error
				conn_.try_recover_from_error();
			}
			catch(shutdown_error &) {
				// termination due to connection shutdown
				throw lost_error("The inlet has been disengaged.");
			}
			catch(std::exception &e) {
				// some perhaps more serious transmission or parsing error (could be indicative of a protocol issue)
				if (!conn_.shutdown())
					std::cerr << "Stream transmission broke off (" << e.what() << "); re-connecting..." << std::endl;
				conn_.try_recover_from_error();
			}
			// wait for a few msec so as to not spam the provider with reconnects
			lslboost::this_thread::sleep_for(lslboost::chrono::milliseconds(500));
		}
	}
	catch(lost_error &) {
		// the connection was irrecoverably lost: since the pull_sample() function may
		// be waiting for the next sample we need to wake it up by passing a sentinel
		sample_queue_.push_sample(sample_p());
	}
	conn_.release_watchdog();
}

} // namespace lsl

max_stars_repo_name: phfix/liblsl

max_issues_repo_name: phfix/liblsl

max_forks_repo_name: phfix/liblsl


text:
#include <string>
#include <pthread.h>
#include <boost/thread/once.hpp>
#include <cstring>

extern "C" int CLogPrintStr(const char *);

int LogPrintStr(const std::string &str) {
  return CLogPrintStr(str.c_str());
}

extern "C"
void plugin_preload_init_cpp() {
  // Trigger global static initialization of the locale system
  std::locale::locale();
}

namespace std  _GLIBCXX_VISIBILITY(default){
  _GLIBCXX_BEGIN_NAMESPACE_VERSION
  void locale::_S_initialize() {
    assert(0);
  }
}

/*
namespace preload {
  extern pthread_key_t &_boosthack_current_thread_tls_key();
  extern boost::once_flag &_boosthack_current_thread_tls_init_flag(void);
  extern pthread_key_t &_boosthack_epoch_tss_key();
  extern pthread_once_t &_boosthack_epoch_tss_key_flag();
  extern BOOST_THREAD_DECL boost::uintmax_t &_boosthack_once_global_epoch();
  extern BOOST_THREAD_DECL pthread_mutex_t &_boosthack_once_epoch_mutex();
  extern BOOST_THREAD_DECL pthread_cond_t &_boosthack_once_epoch_cv();
}

namespace boost {
  namespace detail {
    pthread_key_t &_boosthack_current_thread_tls_key() {
      return preload::_boosthack_current_thread_tls_key();
    }
    boost::once_flag &_boosthack_current_thread_tls_init_flag() {
      return preload::_boosthack_current_thread_tls_init_flag();
    }
    pthread_key_t &_boosthack_epoch_tss_key() {
      return preload::_boosthack_epoch_tss_key();
    }
    pthread_once_t &_boosthack_epoch_tss_key_flag() {
      return preload::_boosthack_epoch_tss_key_flag();
    }
    BOOST_THREAD_DECL boost::uintmax_t &_boosthack_once_global_epoch() {
      return preload::_boosthack_once_global_epoch();
    }
    BOOST_THREAD_DECL pthread_mutex_t &_boosthack_once_epoch_mutex() {
      return preload::_boosthack_once_epoch_mutex();
    }
    BOOST_THREAD_DECL pthread_cond_t &_boosthack_once_epoch_cv() {
      return preload::_boosthack_once_epoch_cv();
    }
  }
}

*/

max_stars_repo_name: amiller/shadow-plugin-bitcoin

max_issues_repo_name: amiller/shadow-plugin-bitcoin

max_forks_repo_name: amiller/shadow-plugin-bitcoin


text:
#include <set>
#include <map>

#include "vtkzCFDReader.h"

#include "vtkAppendCompositeDataLeaves.h"
#include "vtkCharArray.h"
#include "vtkCollection.h"
#include "vtkDataArraySelection.h"
#include "vtkDirectory.h"
#include "vtkDoubleArray.h"
#include "vtkFieldData.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkIntArray.h"
#include "vtkMultiBlockDataSet.h"
#include "vtkUnstructuredGrid.h"
#include "vtkMultiProcessController.h"
#include "vtkObjectFactory.h"
#include "vtkSortDataArray.h"
#include "vtkStdString.h"
#include "vtkStreamingDemandDrivenPipeline.h"
#include "vtkStringArray.h"
#include <vtkCellType.h>
#include <vtkCellLocator.h>
#include <vtkIntArray.h>
#include <vtkFloatArray.h>
#include <vtkDoubleArray.h>
#include <vtkPointData.h>
#include <vtkCellData.h>
#include <vtkPolyData.h>
#include <vtkSmartPointer.h>

#include <vtksys/SystemTools.hxx>

#include <vtk_jsoncpp.h>

#include <hdf5/hdffile.hpp>
#include <hdf5/hdfdataset.hpp>

#include <boost/python.hpp>

namespace detail
{

  struct BoundaryCondition
  {

  enum {
    NONE=0x0,INTERIOR=0x2,WALL=0x3,INFLOW=0x4,OUTFLOW=0x5,SYMMETRY=0x7,FARFIELD=0x9,PERIODIC=0x12
    };

    static std::string toString(int bc);

  };

  std::string BoundaryCondition::toString(int bc)
  {
    std::string name = "UNKNOWN";

    switch(bc)
    {
    case NONE:
    case INTERIOR:
      name = "interior";
      break;
    case WALL:
      name = "wall";
      break;
    case SYMMETRY:
      name = "symmetry";
      break;
    case FARFIELD:
      name = "farfield";
      break;
    case INFLOW:
      name = "inflow";
      break;
    case OUTFLOW:
      name = "outflow";
      break;
    case PERIODIC:
      name = "periodic";
      break;
    }
    return name;
  }

  /** Cell type encoder
   *  The first five bits are used to encode the fundamental cell type
   *  The next 11bits encode the number of faces (max 2047 faces)
   *  The remaining bits encode the order
   */
  struct CellType
  {
    enum{
      TETRA=0x1,PRISM=0x2,PYRA=0x4,HEX=0x8,POLY=0x10,
      NUM_TYPE=0x5,TYPE_MASK=0x1F,NUM_FACE_MASK=0x7FF
    };

    static int getType(int cellType)
    {
      return cellType & TYPE_MASK;
    }

    static int getNumFaces(int cellType)
    {
      return ( cellType >> NUM_TYPE ) & NUM_FACE_MASK;
    }

    static int encode(int cellType, int numFaces)
    {
      return (numFaces << 5) + cellType;
    }

  };

}

vtkStandardNewMacro(vtkzCFDReader);
vtkCxxSetObjectMacro(vtkzCFDReader, Controller, vtkMultiProcessController);

//-----------------------------------------------------------------------------
vtkzCFDReader::vtkzCFDReader()
{

  this->ZoneDataArraySelection = vtkDataArraySelection::New();
  this->CellDataArraySelection = vtkDataArraySelection::New();

  this->FileName = NULL;
  this->CaseFileName = NULL;
  this->FileNameOld = new vtkStdString;
  this->ProblemName = new vtkStdString;
  this->CaseName = new vtkStdString;

  this->Refresh = false;

  this->NumProcesses = 1;
  this->ProcessId = 0;

//#ifdef PARAVIEW_USE_MPI
  this->Controller = NULL;
  this->SetController(vtkMultiProcessController::GetGlobalController());

  if (this->Controller == NULL)
    {
    this->NumProcesses = 1;
    this->ProcessId = 0;
    }
  else
    {
    this->NumProcesses = this->Controller->GetNumberOfProcesses();
    this->ProcessId = this->Controller->GetLocalProcessId();
    }
//#endif
  //this->CaseType = RECONSTRUCTED_CASE;
  this->MTimeOld = 0;

  this->SetNumberOfInputPorts(0);
  this->SetNumberOfOutputPorts(1);

}

//-----------------------------------------------------------------------------
vtkzCFDReader::~vtkzCFDReader()
{
  delete this->FileNameOld;

  delete this->ProblemName;
  delete this->CaseName;

  this->ZoneDataArraySelection->Delete();

  this->SetController(NULL);
}

//-----------------------------------------------------------------------------
void vtkzCFDReader::PrintSelf(ostream &os, vtkIndent indent)
{
  this->Superclass::PrintSelf(os, indent);
  //os << indent << "Case Type: " << this->CaseType << endl;
  os << indent << "MTimeOld: " << this->MTimeOld << endl;
  os << indent << "Number of Processes: " << this->NumProcesses << endl;
  os << indent << "Process Id: " << this->ProcessId << endl;
  os << indent << "Controller: " << this->Controller << endl;
}

//-----------------------------------------------------------------------------
/*
void vtkzCFDReader::SetCaseType(const int t)
{
  if (this->CaseType != t)
    {
    this->CaseType = static_cast<caseType>(t);
    this->Refresh = true;
    this->Modified();
    }
}
*/
//-----------------------------------------------------------------------------

void vtkzCFDReader::ReadPython(std::map<int,std::string> &zoneToBc)
{

  using namespace boost::python;

  const char *env = vtksys::SystemTools::GetEnv("ZCFD_HOME");

  if(env)
  {
    PyObject *obj = PySys_GetObject("path");
    PyObject *pPath = PyUnicode_FromString(env);
    PyList_Append(obj,pPath);
    Py_DECREF(pPath);

    PySys_SetObject("path",obj);

    //PySys_SetPath(const_cast<char*>(env));
  }

  std::cout << "Reading Case File " <<  *CaseName << " " << std::string(env) << std::endl;

  //Initialize python
  Py_Initialize();

  PyRun_SimpleString("import sys\nprint sys.path");

  //Get the main module
  object main_module = import("__main__");
  object main_namespace = main_module.attr("__dict__");

  try
  {
    object perror = exec_file( (const char *)(*CaseName), main_namespace, main_namespace);
  }
  catch(const boost::python::error_already_set&)
  {
    PyObject *e, *v, *t;
    PyErr_Fetch(&e, &v, &t);

    //Get error message
    char *pStrErrorMessage = PyBytes_AsString(v);

    std::cout << "Python Exception: " << pStrErrorMessage << std::endl;

    // A NULL e means that there is not available Python
    // exception
    if (!e) return;
/*
    // See if the exception was an AttributeError. If so,
    // throw a C++ version of that exception
    if (PyErr_GivenExceptionMatches(perror.ptr(), e))
    {
      // We construct objects now since we plan to keep
      // ownership of the references.
      object e_obj(handle<>(allow_null(e));
      object v_obj(handle<>(allow_null(v));
      object t_obj(handle<>(allow_null(t));

      throw AttributeException(e_obj, v_obj, t_obj);
    }
*/
    // We didn't do anything with the Python exception,
    // and we never took ownership of the refs, so it's
    // safe to simply pass them back to PyErr_Restore
    PyErr_Restore(e, v, t);

    // Rethrow the exception (or whatever...this
    // is just an example.)
    throw;
  }

  extract<dict> cppdict_ext(main_namespace["parameters"]);

  assert(cppdict_ext.check());

  dict param = cppdict_ext();

  std::cout << "Read Case File" << std::endl;

  //object r = param.get("reference");
  //std::string s = extract<std::string>(r);
  //std::cout << s;

  bool found = true;
  int count = 1;
  while(found)
  {
    std::ostringstream key;
    key << "BC_" << count;
    if(param.has_key(key.str()))
    {
      object ob = param.get(key.str());
      dict bc = extract<dict>(ob);

      ob = bc.get("type");
      std::string bctype = extract<std::string>(ob);

      if(bc.has_key("zone"))
      {
        object ob = bc.get("zone");
        list l = extract<list>(ob);
        for(int i = 0; i < len(l); ++i)
        {
          object ob = l[i];
          int z = extract<int>(ob);

          zoneToBc[z] = bctype;
        }
      }
    }
    else
      found = false;
    count++;
  }

}

bool vtkzCFDReader::ReadStatusFile(const char *fileName)
{
  std::string status;

  std::string settingsFileName(fileName);
  std::ifstream settingsFile(settingsFileName.c_str(), ios::in | ios::binary | ios::ate);
  vtkDebugMacro("Attempting to load settings file '" << fileName << "'");
  if (settingsFile.is_open())
    {
    std::streampos size = settingsFile.tellg();
    settingsFile.seekg(0, ios::beg);
    int stringSize = size;
    char * settingsString = new char[stringSize+1];
    settingsFile.read(settingsString, stringSize);
    settingsString[stringSize] = '\0';
    settingsFile.close();

    status = settingsString;

    delete[] settingsString;
    }
  std::cout << status << std::endl;

  // Parse the user settings
  Json::Reader reader;
  Json::Value val;
  bool success = reader.parse(status, val, false);

  if(val.isMember("num processor"))
  {
    std::cout << val.get("num processor",Json::Value(0)) << std::endl;
  }
  bool problemNameSet=false;
  if(val.isMember("problem"))
  {

    std::cout << val.get("problem",Json::Value("UNKNOWN")) << std::endl;

    *this->ProblemName = vtksys::SystemTools::GetFilenamePath(fileName) + "/";

    *this->ProblemName += vtkStdString(val.get("problem",Json::Value("UNKNOWN")).asString());

    std::cout << "Problem Name: " << *this->ProblemName << std::endl;

    if(!vtksys::SystemTools::FileExists(*this->ProblemName))
    {
      *this->ProblemName += ".h5";
      problemNameSet = vtksys::SystemTools::FileExists(*this->ProblemName);
    }
    else
      problemNameSet=true;
  }

  if(val.isMember("case"))
  {
    *this->CaseName = vtksys::SystemTools::GetFilenamePath(fileName) + "/";

    *this->CaseName += vtkStdString(val.get("case",Json::Value("UNKNOWN")).asString());

    std::cout << "Case Name: " << *this->CaseName << std::endl;

    if(!vtksys::SystemTools::FileExists(*this->CaseName))
    {
      *this->CaseName += ".py";
      if(!vtksys::SystemTools::FileExists(*this->CaseName))
      {
        *this->CaseName = "";
      }
    }

  }

/*
  for(Json::Value::iterator itr = val.begin(); itr != val.end(); ++itr)
  {
    std::cout << *itr << std::endl;
  }
*/
  return problemNameSet;
}

int vtkzCFDReader::RequestInformation(vtkInformation *vtkNotUsed(request),
    vtkInformationVector **vtkNotUsed(inputVector), vtkInformationVector *outputVector)
{

  const char *env = vtksys::SystemTools::GetEnv("ZCFD_HOME");
  if(env == 0)
  {
    vtkErrorMacro("ZCFD_HOME env has to be specified!");
    return 0;
  }
  std::string zcfdhome(env);
  std::cout << "ZCFD_HOME: " << zcfdhome << std::endl;

  if (!this->FileName || strlen(this->FileName) == 0)
  {
    vtkErrorMacro("FileName has to be specified!");
    return 0;
  }

  if (*FileNameOld != this->FileName || this->Refresh)
  {
    // retain selection status when just refreshing a case
    //if (*this->FileNameOld != "" && *this->FileNameOld != this->FileName)
    {
      // clear selections
      this->ZoneDataArraySelection->RemoveAllArrays();
    }

    std::cout << "RequestInformation: Loading Array selection" << std::endl;

    *this->FileNameOld = vtkStdString(this->FileName);

    // Check if h5 or status file
    std::string extension = vtksys::SystemTools::GetFilenameLastExtension(this->FileName);

    std::cout << "Extension: " << extension << std::endl;

    if(extension == ".txt") // _status.txt
    {
      // Read status file
      bool problemNameSet = ReadStatusFile(this->FileName);
      if(!problemNameSet)
        return 0;
    }
    else if(extension == ".h5") // Mesh file
    {
      *this->ProblemName = vtkStdString(this->FileName);
    }

    vtksys::SystemTools::FileExists(this->FileName);
    vtksys::SystemTools::GetFilenameLastExtension(this->FileName);
    vtksys::SystemTools::GetFilenamePath(FileName);
    vtksys::SystemTools::GetFilenameName(FileName);
    vtksys::SystemTools::GetFilenameWithoutExtension(FileName);

    if(CaseFileName != NULL)
    {
      *this->CaseName = vtkStdString(this->CaseFileName);
      std::cout << "RequestInformation: Case Name " << std::endl;
    }

    //this->InternalMeshSelectionStatus
    //    = this->Parent->GetZonehArrayStatus(internalMeshName.c_str());
    std::map<int,std::string> zoneToBc;
    if(*CaseName != "")
    {
      ReadPython(zoneToBc);
    }

    // Read zone indexes from h5 file
    hdf::HDFFile<> file(*this->ProblemName);
    boost::shared_ptr<hdf::HDFGroup<> > meshg = file.openGroup("mesh",true);
    int totalNumFaces, totalNumCells;

    meshg->openAttribute("numFaces")->readData(totalNumFaces);
    meshg->openAttribute("numCells")->readData(totalNumCells);

    int numCells = totalNumCells;
    int numAllFaces = totalNumFaces;

    std::vector<int> zoneInfo;
    std::vector<int> faceBC;

    faceBC.resize(numAllFaces);
    zoneInfo.resize(numAllFaces*2);

    meshg->openDataset("faceBC")->readData(faceBC);
    meshg->openDataset("faceInfo")->readData(zoneInfo);

    std::map<std::string,std::set<int> > bc;
    std::set<int> zone;
    std::map<int,std::vector<int> > bcToZone;

    assert(faceBC.size() == zoneInfo.size()/2);

    // Count boundary faces
    int numBoundaryFaces = 0;
    for (int i = 0; i < faceBC.size(); i++)
    {
      int z = zoneInfo[i*2];
      if (faceBC[i] != detail::BoundaryCondition::NONE)
      {
        numBoundaryFaces++;
        zone.insert(z);
      }

      if(zoneToBc.find(z) != zoneToBc.end())
      {
        bc[zoneToBc[z]].insert(z);
      }
      else
      {
        std::string str = detail::BoundaryCondition::toString(faceBC[i]);
        if(str != "interior")
          bc[str].insert(z);
      }
    }

    std::set<int> cellZoneSet;
    {
      std::vector<int> cellZones(numCells,0);
      try{
        // Read cells zones

        meshg->openDataset("cellZone")->readData(cellZones);

      }
      catch(...){
      }
      for(int i = 0; i < cellZones.size(); ++i){
        cellZoneSet.insert(cellZones[i]);
      }
    }
    for(std::set<int>::iterator itr = cellZoneSet.begin(); itr != cellZoneSet.end(); ++itr){
      char buf[1024];
      std::sprintf(buf, "%d", *itr);

      bc["interior-"+std::string(buf)].insert(*itr);
    }

    int count = 0;

    selectionToZone.clear();

    std::cout << "RequestInformation: Populating Selection " << GetNumberOfSelectionArrays(ZoneDataArraySelection) << std::endl;

    // For each boundary condition
    for(std::map<std::string,std::set<int> >::iterator it = bc.begin(); it != bc.end(); ++it)
    {
      this->ZoneDataArraySelection->AddArray(it->first.c_str());
      selectionToZone[count] = it->second;
      std::cout << it->first << " " << *(it->second.begin()) << std::endl;
      count++;
    }

    // For each zone create a set
    int numZones = zone.size();
    for(std::set<int>::iterator it = zone.begin(); it != zone.end(); ++it)
    {
      char buf[1024];
      std::sprintf(buf, "%d", *it);
      //this->ZoneDataArraySelection->AddArray(("Zone "+std::to_string(*it)).c_str());
      this->ZoneDataArraySelection->AddArray(("Zone "+std::string(buf)).c_str());
      selectionToZone[count].insert(*it);
      count++;
    }

    double timeRange[2];
    timeRange[0] = timeRange[1] = 0.0;
    outputVector->GetInformationObject(0)->Set(vtkStreamingDemandDrivenPipeline::TIME_STEPS(), timeRange, 0);
    outputVector->GetInformationObject(0)->Set(vtkStreamingDemandDrivenPipeline::TIME_RANGE(), timeRange, 2);

    this->Refresh = false;

  }

  /*

  if (this->CaseType == RECONSTRUCTED_CASE)
    {
    int ret = 1;
    if (this->ProcessId == 0)
      {
      ret = this->Superclass::RequestInformation(request, inputVector,
          outputVector);
      }
    if (this->NumProcesses > 1)
      {
      // if there was an error in process 0 abort all processes
      this->BroadcastStatus(ret);
      if (ret == 0)
        {
        vtkErrorMacro(<< "The master process returned an error.");
        return 0;
        }

      vtkDoubleArray *timeValues;
      if (this->ProcessId == 0)
        {
        timeValues = this->Superclass::GetTimeValues();
        }
      else
        {
        timeValues = vtkDoubleArray::New();
        }
      this->Controller->Broadcast(timeValues, 0);
      if (this->ProcessId != 0)
        {
        this->Superclass::SetTimeInformation(outputVector, timeValues);
        timeValues->Delete();
        this->Superclass::Refresh = false;
        }
      this->GatherMetaData(); // pvserver deadlocks without this
      }

    return ret;
    }

  if (!this->Superclass::FileName || strlen(this->Superclass::FileName) == 0)
    {
    vtkErrorMacro("FileName has to be specified!");
    return 0;
    }

  if (*this->Superclass::FileNameOld != this->Superclass::FileName
      || this->Superclass::ListTimeStepsByControlDict
          != this->Superclass::ListTimeStepsByControlDictOld
      || this->Superclass::Refresh)
    {
    // retain selection status when just refreshing a case
    if (*this->Superclass::FileNameOld != "" && *this->Superclass::FileNameOld != this->Superclass::FileName)
      {
      // clear selections
      this->Superclass::CellDataArraySelection->RemoveAllArrays();
      this->Superclass::PointDataArraySelection->RemoveAllArrays();
      this->Superclass::LagrangianDataArraySelection->RemoveAllArrays();
      this->Superclass::PatchDataArraySelection->RemoveAllArrays();
      }

    *this->Superclass::FileNameOld = vtkStdString(this->FileName);
    this->Superclass::Readers->RemoveAllItems();
    this->Superclass::NumberOfReaders = 0;

    vtkStringArray *procNames = vtkStringArray::New();
    vtkDoubleArray *timeValues;

    // recreate case information
    vtkStdString masterCasePath, controlDictPath;
    this->Superclass::CreateCasePath(masterCasePath, controlDictPath);

    this->Superclass::CreateCharArrayFromString(this->Superclass::CasePath,
        "CasePath", masterCasePath);

    int ret = 1;
    if (this->ProcessId == 0)
      {
      // search and list processor subdirectories
      vtkDirectory *dir = vtkDirectory::New();
      if (!dir->Open(masterCasePath.c_str()))
        {
        vtkErrorMacro(<< "Can't open " << masterCasePath.c_str());
        dir->Delete();
        this->BroadcastStatus(ret = 0);
        return 0;
        }
      vtkIntArray *procNos = vtkIntArray::New();
      for (int fileI = 0; fileI < dir->GetNumberOfFiles(); fileI++)
        {
        const vtkStdString subDir(dir->GetFile(fileI));
        if (subDir.substr(0, 9) == "processor")
          {
          const vtkStdString procNoStr(subDir.substr(9, vtkStdString::npos));
          char *conversionEnd;
          const int procNo = strtol(procNoStr.c_str(), &conversionEnd, 10);
          if (procNoStr.c_str() + procNoStr.length() == conversionEnd && procNo
              >= 0)
            {
            procNos->InsertNextValue(procNo);
            procNames->InsertNextValue(subDir);
            }
          }
        }
      procNos->Squeeze();
      procNames->Squeeze();
      dir->Delete();

      // sort processor subdirectories by processor numbers
      vtkSortDataArray::Sort(procNos, procNames);
      procNos->Delete();

      // get time directories from the first processor subdirectory
      if (procNames->GetNumberOfTuples() > 0)
        {
        vtkOpenFOAMReader *masterReader = vtkOpenFOAMReader::New();
        masterReader->SetFileName(this->FileName);
        masterReader->SetParent(this);
        if (!masterReader->MakeInformationVector(outputVector, procNames
        ->GetValue(0)) || !masterReader->MakeMetaDataAtTimeStep(true))
          {
          procNames->Delete();
          masterReader->Delete();
          this->BroadcastStatus(ret = 0);
          return 0;
          }
        this->Superclass::Readers->AddItem(masterReader);
        timeValues = masterReader->GetTimeValues();
        masterReader->Delete();
        }
      else
        {
        timeValues = vtkDoubleArray::New();
        this->Superclass::SetTimeInformation(outputVector, timeValues);
        }
      }
    else
      {
      timeValues = vtkDoubleArray::New();
      }

    if (this->NumProcesses > 1)
      {
      // if there was an error in process 0 abort all processes
      this->BroadcastStatus(ret);
      if (ret == 0)
        {
        vtkErrorMacro(<< "The master process returned an error.");
        timeValues->Delete(); // don't have to care about process 0
        return 0;
        }

      this->Broadcast(procNames);
      this->Controller->Broadcast(timeValues, 0);
      if (this->ProcessId != 0)
        {
        this->Superclass::SetTimeInformation(outputVector, timeValues);
        timeValues->Delete();
        }
      }

    if (this->ProcessId == 0 && procNames->GetNumberOfTuples() == 0)
      {
      timeValues->Delete();
      }

    // create reader instances for other processor subdirectories
    // skip processor0 since it's already created
    for (int procI = (this->ProcessId ? this->ProcessId : this->NumProcesses); procI
        < procNames->GetNumberOfTuples(); procI += this->NumProcesses)
      {
      vtkOpenFOAMReader *subReader = vtkOpenFOAMReader::New();
      subReader->SetFileName(this->FileName);
      subReader->SetParent(this);
      // if getting metadata failed simply delete the reader instance
      if (subReader->MakeInformationVector(NULL, procNames->GetValue(procI))
          && subReader->MakeMetaDataAtTimeStep(true))
        {
        this->Superclass::Readers->AddItem(subReader);
        }
      else
        {
        vtkWarningMacro(<<"Removing reader for processor subdirectory "
            << procNames->GetValue(procI).c_str());
        }
      subReader->Delete();
      }

    procNames->Delete();

    this->GatherMetaData();
    this->Superclass::Refresh = false;
    }

  outputVector->GetInformationObject(0)->Set(
    CAN_HANDLE_PIECE_REQUEST(),
    1);
*/
  return 1;
}

//-----------------------------------------------------------------------------
int vtkzCFDReader::RequestData(vtkInformation *vtkNotUsed(request),
    vtkInformationVector **vtkNotUsed(inputVector), vtkInformationVector *outputVector)
{

  vtkInformation* outInfo = outputVector->GetInformationObject(0);
  vtkMultiBlockDataSet
      *output =
          vtkMultiBlockDataSet::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT()));

  int nSteps = 0;
  double requestedTimeValue(0);
  if (outInfo->Has(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEP()))
    {
    requestedTimeValue
        = outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEP());
    nSteps = outInfo->Length(vtkStreamingDemandDrivenPipeline::TIME_STEPS());
    }

  // Read request zones
  hdf::HDFFile<> file(*this->ProblemName);
  boost::shared_ptr<hdf::HDFGroup<> > meshg = file.openGroup("mesh",false);
  int totalNumFaces, totalNumCells;

  meshg->openAttribute("numFaces")->readData(totalNumFaces);
  meshg->openAttribute("numCells")->readData(totalNumCells);

  int numCells = totalNumCells;
  int numAllFaces = totalNumFaces;

  std::vector<int> zoneInfo;
  //std::vector<int> faceBC;

  //faceBC.resize(numAllFaces);
  zoneInfo.resize(numAllFaces*2);

  //meshg->openDataset("faceBC")->readData(faceBC);
  meshg->openDataset("faceInfo")->readData(zoneInfo);

  std::vector<int> faceTypeGobalPtr;
  faceTypeGobalPtr.reserve(numAllFaces+1);

  meshg->openDataset("faceType")->readData(faceTypeGobalPtr);
  faceTypeGobalPtr.push_back(0);
  for(int i=faceTypeGobalPtr.size()-1;i>0;--i)
    faceTypeGobalPtr[i] = faceTypeGobalPtr[i-1];
  faceTypeGobalPtr[0] = 0;
  for(int i=1;i<faceTypeGobalPtr.size();++i)
    faceTypeGobalPtr[i] += faceTypeGobalPtr[i-1];

  std::vector<int> faceNodes;
  faceNodes.reserve(faceTypeGobalPtr[numAllFaces]);
  meshg->openDataset("faceNodes")->readData(faceNodes);

  std::vector<double> vertex;
  meshg->openDataset("nodeVertex")->readData(vertex);

  std::vector<int> faceCell(numAllFaces*2);
  meshg->openDataset("faceCell")->readData(faceCell);

  // Read cell zones
  std::vector<int> cellZones(numCells,0);
  try{
    // Read cells zones

    meshg->openDataset("cellZone")->readData(cellZones);

  }
  catch(...){
  }

  std::vector<int> cellType(totalNumCells);
  // Check for cellFace
  try
  {
    meshg->openDataset("cellType")->readData(cellType);
  }
  catch(...)
  {
    std::cout << "zCFDReader: cellType missing. Treating all cells as polyhedral " << std::endl;
    std::vector<hsize_t> dimsf(1);
    dimsf[0] = totalNumCells;
    hdf::Slab<1> filespace(dimsf);
    boost::shared_ptr<hdf::HDFDataSet<> > dataset
      = meshg->createDataset<int>("cellType", filespace);

    //Read cell face
    std::vector<int> cellFaceCount(totalNumCells,0);
    for(int i=0;i<numAllFaces;++i)
    {
      int left = faceCell[i*2];
      int right = faceCell[i*2 + 1];
      if(left < totalNumCells)
        cellFaceCount[left]++;
      if(right < totalNumCells)
        cellFaceCount[right]++;
    }
    for(int i=0;i<totalNumCells;++i)
    {
      cellType[i] = detail::CellType::encode(detail::CellType::POLY,cellFaceCount[i]);
    }
    dataset->writeData(cellType);
  }
  std::vector<int> cellFacePtr(totalNumCells+1,0);
  for(int i=0;i<totalNumCells;++i)
  {
    cellFacePtr[i+1] = cellFacePtr[i] + detail::CellType::getNumFaces(cellType[i]);
  }
  std::vector<int> cellFace(cellFacePtr[totalNumCells]);
  try
  {
    meshg->openDataset("cellFace")->readData(cellFace);
  }
  catch(...)
  {
    // Need to create cellFaces
    std::cout << "zCFDReader: Creating cell face information " << std::endl;
    std::vector<hsize_t> dimsf(1);
    dimsf[0] = cellFace.size();
    hdf::Slab<1> filespace(dimsf);
    boost::shared_ptr<hdf::HDFDataSet<> > dataset
      = meshg->createDataset<int>("cellFace", filespace);

    std::vector<int> count(totalNumCells,0);
    for(int i = 0; i < numAllFaces; ++i)
    {
      for(int j = 0; j < 2; ++j)
      {
        int c = faceCell[i*2 + j];
        if(c < totalNumCells)
        {
          cellFace[cellFacePtr[c]+count[c]] = i;
          count[c]++;
        }
      }
    }
    dataset->writeData(cellFace);
  }

  for(int i=0;i<GetNumberOfSelectionArrays(ZoneDataArraySelection);++i)
  {
    std::string name = GetSelectionArrayName(ZoneDataArraySelection,i);

    if(name.find("interior") != 0 && GetSelectionArrayStatus(ZoneDataArraySelection,name.c_str()))
    {
      // User has selected zone for extraction

      // Get list of faces
      std::set<int> &zoneSet = selectionToZone[i];
      std::vector<int> faceList;
      for(int j=0;j<zoneInfo.size()/2;++j)
      {
        int z = zoneInfo[j*2];
        if(zoneSet.find(z) != zoneSet.end())
        {
          faceList.push_back(j);
          //std::cout << "Zones: " << z << std::endl;
        }
      }
      int numFaces = faceList.size();

      std::cout << "Faces: " << numFaces << std::endl;

      vtkSmartPointer<vtkPolyData> polyData =
          vtkSmartPointer<vtkPolyData>::New();
      polyData->Allocate(numFaces); // set number of faces

      vtkSmartPointer<vtkFloatArray> var =
          vtkSmartPointer<vtkFloatArray>::New();
      var->SetName("zone");
      var->SetNumberOfComponents(1);
      var->SetNumberOfValues(numFaces);

      vtkIdType cellList[100];
      int count = 0;
      int faceCount = 0;
      std::vector<int> nodeMap(vertex.size()/3, -1);
      std::vector<int> reverseMap;

      // Extract faces
      for(int f=0;f<faceList.size();++f)
      {
        int face=faceList[f];
        int z = zoneInfo[face*2];
        int start = faceTypeGobalPtr[face];
        int end   = faceTypeGobalPtr[face+1];
        int numPts = end - start;
        assert(numPts < 100);

        for (int j = start; j < end; ++j)
        {
          int n = faceNodes[j];
          int mapIndex = -1;
          if (nodeMap[n] == -1)
          {
            mapIndex = count;
            nodeMap[n] = count;
            reverseMap.push_back(n);
            count++;
          }
          else
          {
            mapIndex = nodeMap[n];
          }
          cellList[j - start] = mapIndex;
        }

        int vtkCellType = VTK_EMPTY_CELL;
        if (numPts == 3)
          vtkCellType = VTK_TRIANGLE;
        else if (numPts == 4)
          vtkCellType = VTK_QUAD;
        else
          vtkCellType = VTK_POLYGON;

        polyData->InsertNextCell(vtkCellType, numPts, cellList);

        float va = z;
        var->InsertTuple(faceCount,&va);
        faceCount++;
      }

      polyData->GetCellData()->SetActiveScalars("zone");
      polyData->GetCellData()->SetScalars(var);

      vtkSmartPointer<vtkPoints> newPts = vtkSmartPointer<vtkPoints>::New();
      newPts->SetNumberOfPoints(count);

      float xyz[3];
      for (int v = 0; v < 3; ++v)
        xyz[v] = 0.0;

      for (int i = 0; i < reverseMap.size(); ++i)
      {
        int n = reverseMap[i];

        for (int v = 0; v < 3; ++v)
          xyz[v] = vertex[3*n+v];

        newPts->InsertPoint(i, xyz[0], xyz[1], xyz[2]);
      }
      polyData->SetPoints(newPts);

      // Extract nodes
      //vtkUnstructuredGrid *ugrid = vtkUnstructuredGrid::New();

      //vtkMultiBlockDataSet *subOutput = vtkMultiBlockDataSet::New();

      const int blockI = output->GetNumberOfBlocks();
      output->SetBlock(blockI, polyData);
      output->GetMetaData(blockI)->Set(vtkCompositeDataSet::NAME(), name.c_str());
      //polyData->Delete();

    }
    else if(name.find("interior") == 0 && GetSelectionArrayStatus(ZoneDataArraySelection,name.c_str()))
    {

      // Get list of cells
      std::set<int> &zoneSet = selectionToZone[i];

      std::vector<int> cellList;
      cellList.reserve(totalNumCells);
      for(std::set<int>::iterator it = zoneSet.begin(); it != zoneSet.end(); ++it){
        for(int i = 0; i < totalNumCells; ++i){
          if(cellZones[i] == *it){
            cellList.push_back(i);
          }
        }
      }

      // Add interior
      vtkSmartPointer<vtkUnstructuredGrid> uGrid = vtkSmartPointer<
          vtkUnstructuredGrid>::New();
      uGrid->Allocate(cellList.size());

      for(int j = 0; j < cellList.size(); ++j)
      {
        int i = cellList[j];

        std::set<vtkIdType> cellNodeList;
        std::vector<vtkIdType> faceStream;

        int start = cellFacePtr[i];
        int end  = cellFacePtr[i+1];

        for(int f = start; f < end; ++f)
        {
          int face = cellFace[f];
          int z = zoneInfo[face*2];

          int start = faceTypeGobalPtr[face];
          int end   = faceTypeGobalPtr[face+1];
          int numPts = end - start;
          assert(numPts < 100);
          faceStream.push_back(numPts);

          for (int j = start; j < end; ++j)
          {
            int n = faceNodes[j];
            cellNodeList.insert(n);
            faceStream.push_back(n);
          }
        }
        std::vector<vtkIdType> cNodeList(cellNodeList.begin(),
            cellNodeList.end());
        vtkIdType nFaces = end - start;
        vtkIdType nNodes = cNodeList.size();

        uGrid->InsertNextCell(VTK_POLYHEDRON, nNodes, &cNodeList[0], nFaces,
            &faceStream[0]);
      }

      vtkSmartPointer<vtkPoints> newPts = vtkSmartPointer<vtkPoints>::New();
      newPts->SetNumberOfPoints(vertex.size());

      float xyz[3];
      for (int v = 0; v < 3; ++v)
        xyz[v] = 0.0;

      for (int i = 0; i < vertex.size()/3; ++i)
      {
        for (int v = 0; v < 3; ++v)
          xyz[v] = vertex[i*3+v];

        newPts->InsertPoint(i, xyz[0], xyz[1], xyz[2]);
      }
      uGrid->SetPoints(newPts);

      const int blockI = output->GetNumberOfBlocks();
      output->SetBlock(blockI, uGrid);
      output->GetMetaData(blockI)->Set(vtkCompositeDataSet::NAME(), name.c_str());
    }
  }

  /*
  if (this->CaseType == RECONSTRUCTED_CASE)
    {
    int ret = 1;
    if (this->ProcessId == 0)
      {
      ret = this->Superclass::RequestData(request, inputVector, outputVector);
      }
    this->BroadcastStatus(ret);
    this->GatherMetaData();
    return ret;
    }

  vtkInformation* outInfo = outputVector->GetInformationObject(0);
  vtkMultiBlockDataSet
      *output =
          vtkMultiBlockDataSet::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT()));

  int ret = 1;
  if (this->Superclass::Readers->GetNumberOfItems() > 0)
    {
    int nSteps = 0;
    double requestedTimeValue(0);
    if (outInfo->Has(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEP()))
      {
      requestedTimeValue
          = outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEP());
      nSteps = outInfo->Length(vtkStreamingDemandDrivenPipeline::TIME_STEPS());
      if (nSteps > 0)
        {
        outInfo->Set(vtkDataObject::DATA_TIME_STEP(), requestedTimeValue);
        }
      }

    vtkAppendCompositeDataLeaves *append = vtkAppendCompositeDataLeaves::New();
    // append->AppendFieldDataOn();

    vtkOpenFOAMReader *reader;
    this->Superclass::CurrentReaderIndex = 0;
    this->Superclass::Readers->InitTraversal();
    while ((reader
        = vtkOpenFOAMReader::SafeDownCast(this->Superclass::Readers->GetNextItemAsObject()))
        != NULL)
      {
      // even if the child readers themselves are not modified, mark
      // them as modified if "this" has been modified, since they
      // refer to the property of "this"
      if ((nSteps > 0 && reader->SetTimeValue(requestedTimeValue))
          || this->MTimeOld != this->GetMTime())
        {
        reader->Modified();
        }
      if (reader->MakeMetaDataAtTimeStep(false))
        {
        append->AddInputConnection(reader->GetOutputPort());
        }
      }

    this->GatherMetaData();

    if (append->GetNumberOfInputConnections(0) == 0)
      {
      output->Initialize();
      ret = 0;
      }
    else
      {
      // reader->RequestInformation() and RequestData() are called
      // for all reader instances without setting UPDATE_TIME_STEPS
      append->Update();
      output->ShallowCopy(append->GetOutput());
      }
    append->Delete();

    // known issue: output for process without sub-reader will not have CasePath
    output->GetFieldData()->AddArray(this->Superclass::CasePath);
    }
  else
    {
    this->GatherMetaData();
    // page 322 of The ParaView Guide says the output must be initialized
    output->Initialize();
    }

  this->Superclass::UpdateStatus();
  this->MTimeOld = this->GetMTime();
*/
  return 1;
}

int vtkzCFDReader::GetNumberOfSelectionArrays(vtkDataArraySelection *s)
{
  return s->GetNumberOfArrays();
}
int vtkzCFDReader::GetSelectionArrayStatus(vtkDataArraySelection *s,
    const char *name)
{
  return s->ArrayIsEnabled(name);
}
void vtkzCFDReader::SetSelectionArrayStatus(vtkDataArraySelection *s,
    const char* name, int status)
{
  unsigned long int mTime = s->GetMTime();
  if (status)
    {
    s->EnableArray(name);
    }
  else
    {
    s->DisableArray(name);
    }
  if (mTime != s->GetMTime()) // indicate that the pipeline needs to be updated
    {
    this->Modified();
    }
}
const char *vtkzCFDReader::GetSelectionArrayName(vtkDataArraySelection *s,
    int index)
{
  return s->GetArrayName(index);
}
//-----------------------------------------------------------------------------
void vtkzCFDReader::BroadcastStatus(int &status)
{
  if (this->NumProcesses > 1)
    {
    this->Controller->Broadcast(&status, 1, 0);
    }
}

//-----------------------------------------------------------------------------
void vtkzCFDReader::GatherMetaData()
{
  if (this->NumProcesses > 1)
    {
    /*
    this->AllGather(this->Superclass::PatchDataArraySelection);
    this->AllGather(this->Superclass::CellDataArraySelection);
    this->AllGather(this->Superclass::PointDataArraySelection);
    this->AllGather(this->Superclass::LagrangianDataArraySelection);
    // omit removing duplicated entries of LagrangianPaths as well
    // when the number of processes is 1 assuming there's no duplicate
    // entry within a process
    this->AllGather(this->Superclass::LagrangianPaths);
    */
    }
}

//-----------------------------------------------------------------------------
// Broadcast a vtkStringArray in process 0 to all processes
void vtkzCFDReader::Broadcast(vtkStringArray *sa)
{
  vtkIdType lengths[2];
  if (this->ProcessId == 0)
    {
    lengths[0] = sa->GetNumberOfTuples();
    lengths[1] = 0;
    for (int strI = 0; strI < sa->GetNumberOfTuples(); strI++)
      {
      lengths[1] += static_cast<vtkIdType>(sa->GetValue(strI).length()) + 1;
      }
    }
  this->Controller->Broadcast(lengths, 2, 0);
  char *contents = new char [lengths[1]];
  if (this->ProcessId == 0)
    {
    for (int strI = 0, idx = 0; strI < sa->GetNumberOfTuples(); strI++)
      {
      const int len = static_cast<int>(sa->GetValue(strI).length()) + 1;
      memmove(contents + idx, sa->GetValue(strI).c_str(), len);
      idx += len;
      }
    }
  this->Controller->Broadcast(contents, lengths[1], 0);
  if (this->ProcessId != 0)
    {
    sa->Initialize();
    sa->SetNumberOfTuples(lengths[0]);
    for (int strI = 0, idx = 0; strI < lengths[0]; strI++)
      {
      sa->SetValue(strI, contents + idx);
      idx += static_cast<int>(sa->GetValue(strI).length()) + 1;
      }
    }
  delete [] contents;
}

//-----------------------------------------------------------------------------
// AllGather vtkStringArray from and to all processes
void vtkzCFDReader::AllGather(vtkStringArray *s)
{
  vtkIdType length = 0;
  for (int strI = 0; strI < s->GetNumberOfTuples(); strI++)
    {
    length += static_cast<vtkIdType>(s->GetValue(strI).length()) + 1;
    }
  vtkIdType *lengths = new vtkIdType [this->NumProcesses];
  this->Controller->AllGather(&length, lengths, 1);
  vtkIdType totalLength = 0;
  vtkIdType *offsets = new vtkIdType [this->NumProcesses];
  for (int procI = 0; procI < this->NumProcesses; procI++)
    {
    offsets[procI] = totalLength;
    totalLength += lengths[procI];
    }
  char *allContents = new char [totalLength], *contents = new char [length];
  for (int strI = 0, idx = 0; strI < s->GetNumberOfTuples(); strI++)
    {
    const int len = static_cast<int>(s->GetValue(strI).length()) + 1;
    memmove(contents + idx, s->GetValue(strI).c_str(), len);
    idx += len;
    }
  this->Controller->AllGatherV(contents, allContents, length, lengths, offsets);
  delete [] contents;
  delete [] lengths;
  delete [] offsets;
  s->Initialize();
  for (int idx = 0; idx < totalLength; idx += static_cast<int>(strlen(allContents + idx)) + 1)
    {
    const char *str = allContents + idx;
    // insert only when the same string is not found
    if (s->LookupValue(str) == -1)
      {
      s->InsertNextValue(str);
      }
    }
  s->Squeeze();
  delete [] allContents;
}

//-----------------------------------------------------------------------------
// AllGather vtkDataArraySelections from and to all processes
void vtkzCFDReader::AllGather(vtkDataArraySelection *s)
{
  vtkIdType length = 0;
  for (int strI = 0; strI < s->GetNumberOfArrays(); strI++)
    {
    length += static_cast<vtkIdType>(strlen(s->GetArrayName(strI))) + 2;
    }
  vtkIdType *lengths = new vtkIdType [this->NumProcesses];
  this->Controller->AllGather(&length, lengths, 1);
  vtkIdType totalLength = 0;
  vtkIdType *offsets = new vtkIdType [this->NumProcesses];
  for (int procI = 0; procI < this->NumProcesses; procI++)
    {
    offsets[procI] = totalLength;
    totalLength += lengths[procI];
    }
  char *allContents = new char [totalLength], *contents = new char [length];
  for (int strI = 0, idx = 0; strI < s->GetNumberOfArrays(); strI++)
    {
    const char *arrayName = s->GetArrayName(strI);
    contents[idx] = s->ArrayIsEnabled(arrayName);
    const int len = static_cast<int>(strlen(arrayName)) + 1;
    memmove(contents + idx + 1, arrayName, len);
    idx += len + 1;
    }
  this->Controller->AllGatherV(contents, allContents, length, lengths, offsets);
  delete [] contents;
  delete [] lengths;
  delete [] offsets;
  // do not RemoveAllArray so that the previous arrays are preserved
  // s->RemoveAllArrays();
  for (int idx = 0; idx < totalLength; idx += static_cast<int>(strlen(allContents + idx + 1)) + 2)
    {
    const char *arrayName = allContents + idx + 1;
    s->AddArray(arrayName);
    if (allContents[idx] == 0)
      {
      s->DisableArray(arrayName);
      }
    else
      {
      s->EnableArray(arrayName);
      }
    }
  delete [] allContents;
}

ext: cxx

max_stars_repo_name: zCFD/zCFDReader

max_issues_repo_name: zCFD/zCFDReader

max_forks_repo_name: zCFD/zCFDReader