Commit 67b39a07 authored by Kewei Lu's avatar Kewei Lu
Browse files

Revert "add FlowField class for data management(not yet complete)"

This reverts commit 7a709d92.
parent 09351a4e
Pipeline #56144 running with stage
#include "flow_field.h"
FlowField::FlowField() {}
FlowField::~FlowField() {}
void FlowField::ReadData(char* file_name)
{
FILE * pFile;
pFile = fopen(file_name, "rb");
std::cout << "Read the data: " << file_name << std::endl;
if (pFile == NULL) perror ("Error opening file");
else
{
fread(dim_, sizeof(vtkm::Int32), 3, pFile);
std::cout << "Dimension of the data: " << dim_[0] << "," << dim_[1] << "," << dim_[2] << std::endl;
vtkm::Id num_elements = dim_[0] * dim_[1] * dim_[2] * 3;
std::cout << "Number of elements: " << num_elements << std::endl;
float* data = new float[num_elements];
fread(data, sizeof(float), num_elements, pFile);
for (vtkm::Id i = 0; i < num_elements; i++)
{
float x = data[i];
float y = data[++i];
float z = data[++i];
vtkm::Vec<FieldType, 3> vec_data(x, y, z);
this->field.push_back(Normalize(vec_data));
}
//delete data;
}
}
\ No newline at end of file
#ifndef VTKM_FILTERS_FLOWFIELD_H_
#define VTKM_FILTERS_FLOWFIELD_H_
#include <vtkm/cont/DeviceAdapter.h>
#include <vtkm/cont/ArrayHandle.h>
//now that the device adapter is included set a global typedef
//that is the chosen device tag
typedef VTKM_DEFAULT_DEVICE_ADAPTER_TAG DeviceAdapter;
template< typename FieldType>
struct PortalTypes
{
private:
typedef vtkm::cont::ArrayHandle<FieldType> HandleType;
typedef typename HandleType::template ExecutionTypes<DeviceAdapter> ExecutionTypes;
public:
typedef typename ExecutionTypes::Portal Portal;
typedef typename ExecutionTypes::PortalConst PortalConst;
};
/// Vector normalization
template <typename T>
VTKM_EXEC_CONT_EXPORT
T Normalize(T v)
{
float magnitude = sqrt(vtkm::dot(v, v));
if (magnitude == 0)
return vtkm::make_Vec(0.0f, 0.0f, 0.0f);
else
return 1.0f / magnitude * v;
}
template<typename FieldType>
class FlowField {
public:
FlowField() {}
~FlowField() {}
void ReadData(char* file_name)
{
FILE * pFile;
pFile = fopen(file_name, "rb");
std::cout << "Read the data: " << file_name << std::endl;
if (pFile == NULL) perror ("Error opening file");
else
{
fread(dim_, sizeof(vtkm::Int32), 3, pFile);
std::cout << "Dimension of the data: " << dim_[0] << "," << dim_[1] << "," << dim_[2] << std::endl;
vtkm::Id num_elements = dim_[0] * dim_[1] * dim_[2] * 3;
std::cout << "Number of elements: " << num_elements << std::endl;
float* data = new float[num_elements];
fread(data, sizeof(float), num_elements, pFile);
for (vtkm::Id i = 0; i < num_elements; i++)
{
float x = data[i];
float y = data[++i];
float z = data[++i];
vtkm::Vec<FieldType, 3> vec_data(x, y, z);
this->field_.push_back(Normalize(vec_data));
}
}
std::cout << "size of field: " << this->field_.size() << std::endl;
this->field_array_ = vtkm::cont::make_ArrayHandle(&field_[0], field_.size());
field_array_protal_const_ = field_array_.GetPortalConstControl();
}
//get the velocity at index(x,y,z)
//
VTKM_EXEC_EXPORT
vtkm::Vec<FieldType, 3> GetVel(vtkm::Id3 index, const vtkm::Id &xdim, const vtkm::Id &ydim, const vtkm::Id &zdim)
{
vtkm::Id idx = index[2] * ydim * xdim + index[1] * ydim + index[0];
return field_array_protal_const_.Get(idx);
}
//get the velocity at position(x,y,z), need trilinear interpolation
//
VTKM_EXEC_EXPORT
vtkm::Vec<FieldType, 3> AtPhys(vtkm::Vec<vtkm::Float32, 3> pos, const vtkm::Id &xdim, const vtkm::Id &ydim, const vtkm::Id &zdim)
{
//get eight corner index
vtkm::Id3 idx000, idx001, idx010, idx011, idx100, idx101, idx110, idx111;
idx000[0] = vtkm::Id(floor(pos[0]));
idx000[1] = vtkm::Id(floor(pos[1]));
idx000[2] = vtkm::Id(floor(pos[2]));
idx001 = idx000; idx001[0] = (idx001[0] + 1) <= xdim - 1 ? idx001[0] + 1 : xdim - 1;
idx010 = idx000; idx010[1] = (idx010[1] + 1) <= ydim - 1 ? idx010[1] + 1 : ydim - 1;
idx011 = idx010; idx011[0] = (idx011[0] + 1) <= xdim - 1 ? idx011[0] + 1 : xdim - 1;
idx100 = idx000; idx100[2] = (idx100[2] + 1) <= zdim - 1 ? idx100[2] + 1 : zdim - 1;
idx101 = idx100; idx101[0] = (idx101[0] + 1) <= xdim - 1 ? idx101[0] + 1 : xdim - 1;
idx110 = idx100; idx110[1] = (idx110[1] + 1) <= ydim - 1 ? idx110[1] + 1 : ydim - 1;
idx111 = idx110; idx111[0] = (idx111[0] + 1) <= xdim - 1 ? idx111[0] + 1 : xdim - 1;
//get velocity
vtkm::Vec<FieldType, 3> v000, v001, v010, v011, v100, v101, v110, v111;
v000 = GetVel(idx000, xdim, ydim, zdim);
v001 = GetVel(idx001, xdim, ydim, zdim);
v010 = GetVel(idx010, xdim, ydim, zdim);
v011 = GetVel(idx011, xdim, ydim, zdim);
v100 = GetVel(idx100, xdim, ydim, zdim);
v101 = GetVel(idx101, xdim, ydim, zdim);
v110 = GetVel(idx110, xdim, ydim, zdim);
v111 = GetVel(idx111, xdim, ydim, zdim);
//interpolation
vtkm::Vec<FieldType, 3> v00, v01, v10, v11;
float a = pos[0] - floor(pos[0]);
v00[0] = (1.0f - a) * v000[0] + a * v001[0]; v00[1] = (1.0f - a) * v000[1] + a * v001[1]; v00[2] = (1.0f - a) * v000[2] + a * v001[2];
v01[0] = (1.0f - a) * v010[0] + a * v011[0]; v01[1] = (1.0f - a) * v010[1] + a * v011[1]; v01[2] = (1.0f - a) * v010[2] + a * v011[2];
v10[0] = (1.0f - a) * v100[0] + a * v101[0]; v10[1] = (1.0f - a) * v100[1] + a * v101[1]; v10[2] = (1.0f - a) * v100[2] + a * v101[2];
v11[0] = (1.0f - a) * v110[0] + a * v111[0]; v11[1] = (1.0f - a) * v110[1] + a * v111[1]; v11[2] = (1.0f - a) * v110[2] + a * v111[2];
vtkm::Vec<FieldType, 3> v0, v1;
a = pos[1] - floor(pos[1]);
v0[0] = (1.0f - a) * v00[0] + a * v01[0]; v0[1] = (1.0f - a) * v00[1] + a * v01[1]; v0[2] = (1.0f - a) * v00[2] + a * v01[2];
v1[0] = (1.0f - a) * v10[0] + a * v11[0]; v1[1] = (1.0f - a) * v10[1] + a * v11[1]; v1[2] = (1.0f - a) * v10[2] + a * v11[2];
vtkm::Vec<FieldType, 3> v;
a = pos[2] - floor(pos[2]);
v[0] = (1.0f - a) * v0[0] + v1[0]; v[1] = (1.0f - a) * v0[1] + v1[1]; v[2] = (1.0f - a) * v0[2] + v1[2];
return v;
}
vtkm::Int32 GetXDim() {return dim_[0];};
vtkm::Int32 GetYDim() {return dim_[1];};
vtkm::Int32 GetZDim() {return dim_[2];};
std::vector<vtkm::Vec<FieldType, 3> > GetField() {return this->field_;}
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3> > GetFieldArrayHandle() {return this->field_array_;}
private:
vtkm::Int32 dim_[3];
std::vector<vtkm::Vec<FieldType, 3> > field_;
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3> > field_array_;
typedef typename PortalTypes< vtkm::Vec<FieldType, 3> >::PortalConst FieldPortalType;
FieldPortalType field_array_protal_const_;
};
#endif //VTKM_FILTERS_FLOWFIELD_H_
\ No newline at end of file
Supports Markdown
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