Skip to content
Snippets Groups Projects
Commit 0a529a0f authored by Ali Can Demiralp's avatar Ali Can Demiralp
Browse files

Progress.

parent b06274bd
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
include/dpa/types/regular_grid.hpp
+ 171
− 145
View file @ 0a529a0f
......@@ -19,37 +19,38 @@ namespace dpa
template <typename _element_type, std::size_t _dimensions>
struct regular_grid
{
static constexpr auto dimensions = _dimensions;
static constexpr auto dimensions = _dimensions;
using element_type = _element_type;
using element_type = _element_type;
using complex_type = ???
using domain_type = typename vector_traits<scalar, dimensions>::type;
using index_type = std::array<std::size_t, dimensions>;
using container_type = boost::multi_array<element_type, dimensions>;
using array_view_type = typename container_type::template array_view<dimensions>::type;
using type = regular_grid<element_type, dimensions>;
using gradient_type = regular_grid<typename gradient_traits <element_type, dimensions>::type, dimensions>;
using potential_type = regular_grid<typename potential_traits<element_type, dimensions>::type, dimensions>;
using second_gradient_type = regular_grid<typename gradient_traits <typename gradient_traits <element_type, dimensions>::type, dimensions>::type, dimensions>;
using second_potential_type = regular_grid<typename potential_traits<typename potential_traits<element_type, dimensions>::type, dimensions>::type, dimensions>;
using eigen_decomposition_type = regular_grid<std::pair<element_type, typename potential_traits<element_type, dimensions>::type>, dimensions>;
// Ducks [] on the domain_type.
index_type cell_index (const domain_type& position) const
using domain_type = typename vector_traits<scalar, dimensions>::type;
using index_type = std::array<std::size_t, dimensions>;
using container_type = boost::multi_array<element_type, dimensions>;
using array_view_type = typename container_type::template array_view<dimensions>::type;
using type = regular_grid<element_type, dimensions>;
using gradient_type = regular_grid<typename gradient_traits <element_type, dimensions>::type, dimensions>;
using potential_type = regular_grid<typename potential_traits<element_type, dimensions>::type, dimensions>;
using second_gradient_type = regular_grid<typename gradient_traits <typename gradient_traits <element_type, dimensions>::type, dimensions>::type, dimensions>;
using second_potential_type = regular_grid<typename potential_traits<typename potential_traits<element_type, dimensions>::type, dimensions>::type, dimensions>;
???using eigen_decomposition_type = std::pair<regular_grid<complex_type, dimensions>, regular_grid<typename potential_traits<element_type, dimensions>::complex_type, dimensions>>;
???using symmetric_eigen_decomposition_type = std::pair<regular_grid<element_type, dimensions>, regular_grid<typename potential_traits<element_type, dimensions>::type, dimensions>>;
/// Access.
index_type cell_index (const domain_type& position) const
{
index_type index;
for (std::size_t i = 0; i < dimensions; ++i)
index[i] = std::floor((position[i] - offset[i]) / spacing[i]);
return index;
}
// Ducks [] on the domain_type.
element_type& cell (const domain_type& position)
element_type& cell (const domain_type& position)
{
return data(cell_index(position));
}
// Ducks [] on the domain_type.
bool contains (const domain_type& position) const
bool contains (const domain_type& position) const
{
for (std::size_t i = 0; i < dimensions; ++i)
{
......@@ -59,8 +60,7 @@ struct regular_grid
}
return true;
}
// Ducks [] on the domain_type.
element_type interpolate (const domain_type& position) const
element_type interpolate (const domain_type& position) const
{
domain_type weights ;
index_type start_index;
......@@ -84,8 +84,10 @@ struct regular_grid
intermediates[j] = (scalar(1) - weights[i]) * intermediates[2 * j] + weights[i] * intermediates[2 * j + 1];
return intermediates[0];
}
void apply (std::function<void(const index_type&, element_type&)> function)
/// Iteration.
void apply (std::function<void(const index_type&, element_type&)> function)
{
index_type start_index; start_index.fill(0);
index_type end_index ;
......@@ -94,7 +96,7 @@ struct regular_grid
end_index[i] = data.shape()[i];
permute_for<index_type>([&] (const index_type& index) { function(index, data(index)); }, start_index, end_index, increment);
}
void apply_parallel (std::function<void(const index_type&, element_type&)> function)
void apply_parallel (std::function<void(const index_type&, element_type&)> function)
{
index_type start_index; start_index.fill(0);
index_type end_index ;
......@@ -103,21 +105,21 @@ struct regular_grid
end_index[i] = data.shape()[i];
parallel_permute_for<index_type>([&] (const index_type& index) { function(index, data(index)); }, start_index, end_index, increment);
}
void apply_window (std::function<void(const index_type&, element_type&, array_view_type&)> function, const index_type& window_size)
void apply_window (std::function<void(const index_type&, element_type&, array_view_type&)> function, const index_type& window_size)
{
apply([&] (const index_type& index, element_type& element)
{
apply_window_internal(function, window_size, index, element);
});
}
void apply_window_parallel (std::function<void(const index_type&, element_type&, array_view_type&)> function, const index_type& window_size)
void apply_window_parallel (std::function<void(const index_type&, element_type&, array_view_type&)> function, const index_type& window_size)
{
apply_parallel([&] (const index_type& index, element_type& element)
{
apply_window_internal(function, window_size, index, element);
});
}
void apply_window_internal (std::function<void(const index_type&, element_type&, array_view_type&)> function, const index_type& window_size, const index_type& index, element_type& element)
void apply_window_internal (std::function<void(const index_type&, element_type&, array_view_type&)> function, const index_type& window_size, const index_type& index, element_type& element)
{
boost::detail::multi_array::index_gen<dimensions, dimensions> indices;
for (auto dimension = 0; dimension < dimensions; ++dimension)
......@@ -137,7 +139,10 @@ struct regular_grid
function(index, element, data[indices]);
}
gradient_type gradient (const bool normalize = true)
/// Derivatives/integrals through central differences.
// Scalar/vector operator.
gradient_type gradient (const bool normalize = true)
{
auto& shape = reinterpret_cast<index_type const&>(*data.shape());
auto two_spacing = domain_type(2 * spacing);
......@@ -161,7 +166,13 @@ struct regular_grid
return gradient;
}
potential_type potential () const
// Scalar operator.
second_gradient_type second_gradient (const bool normalize = true)
{
return gradient(normalize).gradient(normalize);
}
// Vector/tensor operator.
potential_type potential () const
{
auto& shape = reinterpret_cast<index_type const&>(*data.shape());
auto two_spacing = domain_type(2 * spacing);
......@@ -177,7 +188,6 @@ struct regular_grid
{
auto partial_start_index = start_index; partial_start_index[dimension] = serial_index ;
auto partial_end_index = end_index ; partial_end_index [dimension] = serial_index + 1;
parallel_permute_for<index_type>([&] (const index_type& index)
{
auto prev_index = index, next_index = index;
......@@ -195,41 +205,28 @@ struct regular_grid
}
return potential;
}
// Scalar-only operator.
second_gradient_type hessian ()
// Tensor operator.
second_potential_type second_potential () const
{
return gradient().gradient(); // In generalized formulation of the gradient, transpose is omitted.
return potential().potential();
}
// Scalar-only operator.
type laplacian ()
{
type laplacian {type::container_type(reinterpret_cast<index_type const&>(*data.shape())), offset, size, spacing};
auto hessian_grid = hessian();
hessian_grid.apply_parallel([&] (const index_type& index, typename second_gradient_type::element_type& element)
{
laplacian.data(index) = element.trace();
});
return laplacian;
}
// Scalar-only operator.
second_gradient_type structure_tensor (const index_type& window_size)
/// Moments.
// Vector operator.
???gradient_type structure_tensor (const index_type& window_size)
{
auto gradient_grid = gradient();
second_gradient_type outer_product_grid {second_gradient_type::container_type(reinterpret_cast<index_type const&>(*data.shape())), offset, size, spacing};
gradient_grid .apply_parallel ([&] (const index_type& index, typename gradient_type::element_type& element)
gradient_type outer_product_grid {gradient_type::container_type(reinterpret_cast<index_type const&>(*data.shape())), offset, size, spacing};
apply_parallel([&] (const index_type& index, typename gradient_type::element_type& element)
{
outer_product_grid.data(index) = element.transpose().eval() * element;
});
second_gradient_type structure_tensor {second_gradient_type::container_type(reinterpret_cast<index_type const&>(*data.shape())), offset, size, spacing};
outer_product_grid.apply_window_parallel([&] (const index_type& index, typename second_gradient_type::element_type& element, typename second_gradient_type::array_view_type& elements)
gradient_type structure_tensor {gradient_type::container_type(reinterpret_cast<index_type const&>(*data.shape())), offset, size, spacing};
outer_product_grid.apply_window_parallel([&] (const index_type& index, typename gradient_type::element_type& element, typename gradient_type::array_view_type& elements)
{
structure_tensor.data(index).setZero();
dpa::iterate(elements, [&] (const typename second_gradient_type::element_type& iteratee)
dpa::iterate(elements, [&] (const typename gradient_type::element_type& iteratee)
{
structure_tensor.data(index).array() += (iteratee.array() / elements.num_elements()); // TODO: Adjustable weights instead of 1/elements?
});
......@@ -237,75 +234,55 @@ struct regular_grid
return structure_tensor;
}
// Vector-only operator.
potential_type divergence ()
{
potential_type divergence {potential_type::container_type(reinterpret_cast<index_type const&>(*data.shape())), offset, size, spacing};
auto gradient_grid = gradient();
gradient_grid.apply_parallel([&] (const index_type& index, typename gradient_type::element_type& element)
{
divergence.data(index) = element.trace();
});
return divergence;
}
// Vector-only operator.
type vorticity ()
/// Convenience.
// Vector/tensor operator.
void normalize ()
{
type vorticity {type::container_type(reinterpret_cast<index_type const&>(*data.shape())), offset, size, spacing};
auto gradient_grid = gradient();
gradient_grid.apply_parallel([&] (const index_type& index, typename gradient_type::element_type& element)
apply_parallel([ ] (const index_type& index, element_type& element)
{
vorticity.data(index)[0] = element(2, 1) - element(1, 2);
vorticity.data(index)[1] = element(0, 2) - element(2, 0);
vorticity.data(index)[2] = element(1, 0) - element(0, 1);
scalar norm = element.norm();
if (norm > std::numeric_limits<scalar>::epsilon())
element /= norm;
});
return vorticity;
}
// Vector-only operator.
potential_type q_criterion ()
// Tensor operator.
second_potential_type trace ()
{
potential_type q_criterion {potential_type::container_type(reinterpret_cast<index_type const&>(*data.shape())), offset, size, spacing};
auto gradient_grid = gradient();
gradient_grid.apply_parallel([&] (const index_type& index, typename gradient_type::element_type& element)
second_potential_type trace {second_potential_type::container_type(reinterpret_cast<index_type const&>(*data.shape())), offset, size, spacing};
apply_parallel([&] (const index_type& index, element_type& element)
{
q_criterion.data(index) =
- (element(0, 0) * element(0, 0) + element(1, 1) * element(1, 1) + element(2, 2)) / scalar(2)
- (element(0, 1) * element(1, 0) + element(0, 2) * element(2, 0) + element(1, 2) * element(2, 1));
trace.data(index) = element.trace();
});
return q_criterion;
return trace;
}
// Vector/tensor-only operator.
void normalize ()
// Tensor operator.
??eigen_decomposition_type eigen_decomposition ()
{
apply_parallel([ ] (const index_type& index, element_type& element)
eigen_decomposition_type eigen_decomposition {eigen_decomposition_type::container_type(reinterpret_cast<index_type const&>(*data.shape())), offset, size, spacing};
apply_parallel([&] (const index_type& index, element_type& element)
{
scalar norm = element.norm();
if (norm > std::numeric_limits<scalar>::epsilon())
element /= norm;
Eigen::EigenSolver<element_type> solver(element);
eigen_decomposition.data(index) = std::make_pair(solver.eigenvectors(), solver.eigenvalues());
});
return eigen_decomposition;
}
// Tensor-only operator.
eigen_decomposition_type eigen_decomposition ()
// Tensor operator.
??symmetric_eigen_decomposition_type symmetric_eigen_decomposition ()
{
eigen_decomposition_type eigen_decomposition {eigen_decomposition_type::container_type(reinterpret_cast<index_type const&>(*data.shape())), offset, size, spacing};
apply_parallel([&] (const index_type& index, element_type& element)
{
Eigen::EigenSolver<element_type> solver(element);
eigen_decomposition.data(index) = std::make_pair(solver.eigenvectors().real(), solver.eigenvalues().real());
// TODO: Sort descending.
Eigen::SelfAdjointEigenSolver<element_type> solver(element);
eigen_decomposition.data(index) = std::make_pair(solver.eigenvectors(), solver.eigenvalues());
});
return eigen_decomposition;
}
// Tensor-only operator.
eigen_decomposition_type reoriented_eigen_decomposition(const index_type& window_size)
// Tensor operator. Requires symmetry.
??symmetric_eigen_decomposition_type reoriented_eigen_decomposition(const index_type& window_size)
{
auto grid = eigen_decomposition();
auto grid = symmetric_eigen_decomposition();
grid.apply_window([&] (const index_type& index, typename eigen_decomposition_type::element_type& element, typename eigen_decomposition_type::array_view_type& elements)
{
auto counter = 0;
......@@ -329,76 +306,125 @@ struct regular_grid
}, window_size);
return grid;
}
// Tensor-only operator.
second_potential_type fractional_anisotropy ()
/// Special operators.
// Scalar operator.
second_gradient_type hessian ()
{
return second_gradient(); // In generalized formulation of the gradient, transpose is omitted.
}
// Scalar operator.
??type laplacian ()
{
return second_gradient().trace();
}
// Tensor operator. Applied to a vector field gradient.
??second_potential_type divergence ()
{
return trace();
}
// Tensor operator. Applied to a 3x3 vector field gradient.
??potential_type vorticity ()
{
potential_type vorticity {potential_type::container_type(reinterpret_cast<index_type const&>(*data.shape())), offset, size, spacing};
apply_parallel([&] (const index_type& index, element_type& element)
{
vorticity.data(index)[0] = element(2, 1) - element(1, 2);
vorticity.data(index)[1] = element(0, 2) - element(2, 0);
vorticity.data(index)[2] = element(1, 0) - element(0, 1);
});
return vorticity;
}
// Tensor operator. Applied to a 3x3 vector field gradient.
??potential_type q_criterion ()
{
potential_type q_criterion {potential_type::container_type(reinterpret_cast<index_type const&>(*data.shape())), offset, size, spacing};
auto gradient_grid = gradient();
gradient_grid.apply_parallel([&] (const index_type& index, typename gradient_type::element_type& element)
{
q_criterion.data(index) =
- (element(0, 0) * element(0, 0) + element(1, 1) * element(1, 1) + element(2, 2)) / scalar(2)
- (element(0, 1) * element(1, 0) + element(0, 2) * element(2, 0) + element(1, 2) * element(2, 1));
});
return q_criterion;
}
??void vector_laplacian()
{
}
// Tensor operator. Requires symmetry.
second_potential_type axial_diffusivity ()
{
second_potential_type output {second_potential_type::container_type(reinterpret_cast<index_type const&>(*data.shape())), offset, size, spacing};
apply_parallel([&] (const index_type& index, element_type& element)
{
auto lambda = element.eigenvalues();
auto mean = lambda.sum() / lambda.size();
output.data(index) = std::sqrt(3.0 / 2.0) * (lambda - mean).norm() / lambda.norm();
auto lambda = element.eigenvalues().real();
std::sort(lambda.data(), lambda.data() + lambda.size(), std::greater<typename second_potential_type::element_type>());
output.data(index) = lambda[0];
});
return output;
}
// Tensor-only operator.
second_potential_type relative_anisotropy ()
// Tensor operator. Requires symmetry.
second_potential_type radial_diffusivity ()
{
second_potential_type output {second_potential_type::container_type(reinterpret_cast<index_type const&>(*data.shape())), offset, size, spacing};
apply_parallel([&] (const index_type& index, element_type& element)
{
auto lambda = element.eigenvalues();
auto mean = lambda.sum() / lambda.size();
output.data(index) = (lambda - mean).norm() / std::sqrt(3) * mean;
auto lambda = element.eigenvalues().real();
std::sort(lambda.data(), lambda.data() + lambda.size(), std::greater<typename second_potential_type::element_type>());
output.data(index) = (lambda.sum() - lambda[0]) / 2;
});
return output;
}
// Tensor-only operator.
second_potential_type volume_ratio ()
// Tensor operator. Requires symmetry.
second_potential_type mean_diffusivity ()
{
second_potential_type output {second_potential_type::container_type(reinterpret_cast<index_type const&>(*data.shape())), offset, size, spacing};
apply_parallel([&] (const index_type& index, element_type& element)
{
auto lambda = element.eigenvalues().real();
output.data(index) = lambda.sum() / lambda.size();
});
return output;
}
// Tensor-only operator.
second_potential_type axial_diffusivity ()
// Tensor operator. Requires symmetry.
second_potential_type fractional_anisotropy ()
{
second_potential_type axial_diffusivity {second_potential_type::container_type(reinterpret_cast<index_type const&>(*data.shape())), offset, size, spacing};
second_potential_type output {second_potential_type::container_type(reinterpret_cast<index_type const&>(*data.shape())), offset, size, spacing};
apply_parallel([&] (const index_type& index, element_type& element)
{
auto lambda = element.eigenvalues();
axial_diffusivity.data(index) =
std::sqrt(0.5)
* std::sqrt(std::pow(lambda[0] - lambda[1], 2) + std::pow(lambda[1] - lambda[2], 2) + std::pow(lambda[2] - lambda[0], 2))
/ std::sqrt(std::pow(lambda[0] , 2) + std::pow(lambda[1] , 2) + std::pow(lambda[2] , 2));
auto lambda = element.eigenvalues().real();
auto mean = lambda.sum() / lambda.size();
output.data(index) = std::sqrt(lambda.size() / scalar(2)) * potential_type::element_type(lambda.array() - mean).norm() / lambda.norm();
});
return axial_diffusivity;
return output;
}
// Tensor-only operator.
second_potential_type mean_diffusivity ()
// Tensor operator. Requires symmetry.
second_potential_type relative_anisotropy ()
{
second_potential_type mean_diffusivity {second_potential_type::container_type(reinterpret_cast<index_type const&>(*data.shape())), offset, size, spacing};
second_potential_type output {second_potential_type::container_type(reinterpret_cast<index_type const&>(*data.shape())), offset, size, spacing};
apply_parallel([&] (const index_type& index, element_type& element)
{
auto lambda = element.eigenvalues();
mean_diffusivity.data(index) =
std::sqrt(0.5)
* std::sqrt(std::pow(lambda[0] - lambda[1], 2) + std::pow(lambda[1] - lambda[2], 2) + std::pow(lambda[2] - lambda[0], 2))
/ std::sqrt(std::pow(lambda[0] , 2) + std::pow(lambda[1] , 2) + std::pow(lambda[2] , 2));
auto lambda = element.eigenvalues().real();
auto mean = lambda.sum() / lambda.size();
output.data(index) = potential_type::element_type(lambda.array() - mean).norm() / (std::sqrt(lambda.size()) * mean);
});
return mean_diffusivity;
return output;
}
// Tensor-only operator.
second_potential_type relative_diffusivity ()
// Tensor operator. Requires symmetry.
second_potential_type volume_ratio ()
{
second_potential_type relative_diffusivity {second_potential_type::container_type(reinterpret_cast<index_type const&>(*data.shape())), offset, size, spacing};
second_potential_type output {second_potential_type::container_type(reinterpret_cast<index_type const&>(*data.shape())), offset, size, spacing};
apply_parallel([&] (const index_type& index, element_type& element)
{
auto lambda = element.eigenvalues();
relative_diffusivity.data(index) =
std::sqrt(0.5)
* std::sqrt(std::pow(lambda[0] - lambda[1], 2) + std::pow(lambda[1] - lambda[2], 2) + std::pow(lambda[2] - lambda[0], 2))
/ std::sqrt(std::pow(lambda[0] , 2) + std::pow(lambda[1] , 2) + std::pow(lambda[2] , 2));
auto lambda = element.eigenvalues().real();
auto mean = lambda.sum() / lambda.size();
output.data(index) = lambda.prod() / std::pow(mean, lambda.size());
});
return relative_diffusivity;
return output;
}
container_type data {};
......
source/main.cpp
+ 51
− 0
View file @ 0a529a0f
#include <cstdint>
#include <iostream>
#include <dpa/pipeline.hpp>
#include <dpa/types/regular_fields.hpp>
std::int32_t main(std::int32_t argc, char** argv)
{
auto scalars = dpa::regular_scalar_field_2d {boost::multi_array<float, 2>(boost::extents[3][3]), dpa::vector2(0, 0), dpa::vector2(1, 1), dpa::vector2(0.5, 0.5)};
scalars.apply_window_parallel([ ] (const std::array<std::size_t, 2>& index, dpa::scalar& element, dpa::regular_scalar_field_2d::array_view_type& elements)
{
element = elements.num_elements();
}, std::array<std::size_t, 2>{3, 3});
auto laplacian = scalars.laplacian ();
auto structure_tensor = scalars.structure_tensor (std::array<std::size_t, 2>{3, 3});
auto vectors = scalars.gradient ();
auto divergence = vectors.divergence ();
auto vorticity = vectors.vorticity ();
auto q_criterion = vectors.q_criterion ();
auto tensors = vectors.gradient ();
auto decomposition = tensors.eigen_decomposition ();
auto reoriented_decomposition = tensors.reoriented_eigen_decomposition(std::array<std::size_t, 2>{3, 3});
auto fa = tensors.fractional_anisotropy ();
auto ra = tensors.relative_anisotropy ();
auto vr = tensors.volume_ratio ();
auto ad = tensors.axial_diffusivity ();
auto md = tensors.mean_diffusivity ();
auto rd = tensors.relative_diffusivity ();
tensors.apply_parallel([&] (const std::array<std::size_t, 2>& index, dpa::matrix2& element)
{
for (auto i = 0; i < 2; ++i)
element.col(i) = reoriented_decomposition.data(index).first.col(i) * reoriented_decomposition.data(index).second[i];
});
auto shape = scalars.data.shape();
for (auto y = 0; y < shape[1]; ++y)
{
for (auto x = 0; x < shape[0]; ++x)
std::cout << scalars.data[x][y] << " ";
std::cout << "\n";
}
for (auto y = 0; y < shape[1]; ++y)
{
for (auto x = 0; x < shape[0]; ++x)
std::cout << vectors.data[x][y] << " ";
std::cout << "\n";
}
for (auto y = 0; y < shape[1]; ++y)
{
for (auto x = 0; x < shape[0]; ++x)
std::cout << tensors.data[x][y] << "\n";
std::cout << "\n";
}
return dpa::pipeline().run(argc, argv);
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment