mirror of
https://gitlab.kitware.com/vtk/vtk-m
synced 2024-09-08 21:33:55 +00:00
Adding MIR filter
-- Adding MIR filter and tables for material reconstruction
This commit is contained in:
Abhishek Yenpure
parent
82f5b5ebcf
commit
7f4b66d057
@ -86,6 +86,7 @@ set(extra_headers
|
||||
Mask.h
|
||||
MaskPoints.h
|
||||
MeshQuality.h
|
||||
MIRFilter.h
|
||||
NDEntropy.h
|
||||
NDHistogram.h
|
||||
ParticleDensityBase.h
|
||||
@ -141,6 +142,7 @@ set(extra_header_template_sources
|
||||
Mask.hxx
|
||||
MaskPoints.hxx
|
||||
MeshQuality.hxx
|
||||
MIRFilter.hxx
|
||||
NDEntropy.hxx
|
||||
NDHistogram.hxx
|
||||
ParticleDensityCloudInCell.hxx
|
||||
|
||||
147
vtkm/filter/MIRFilter.h
Normal file
147
vtkm/filter/MIRFilter.h
Normal file
@ -0,0 +1,147 @@
|
||||
//============================================================================
|
||||
// Copyright (c) Kitware, Inc.
|
||||
// All rights reserved.
|
||||
// See LICENSE.txt for details.
|
||||
//
|
||||
// This software is distributed WITHOUT ANY WARRANTY; without even
|
||||
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
|
||||
// PURPOSE. See the above copyright notice for more information.
|
||||
//============================================================================
|
||||
|
||||
#ifndef vtkm_m_filter_MIRFilter_h
|
||||
#define vtkm_m_filter_MIRFilter_h
|
||||
|
||||
#include <vtkm/cont/ArrayCopy.h>
|
||||
|
||||
#include <vtkm/filter/FilterDataSet.h>
|
||||
#include <vtkm/filter/FilterDataSetWithField.h>
|
||||
|
||||
namespace vtkm
|
||||
{
|
||||
namespace filter
|
||||
{
|
||||
/// @brief Calculates and subdivides a mesh based on the material interface reconstruction algorithm.
|
||||
///
|
||||
/// Subdivides a mesh given volume fraction information for each _cell_. It does this by applying a
|
||||
/// mixture of the painters algorithm and isosurfacing algorithm. This filter will return
|
||||
/// a dataset where cells are subdivided into new cells of a certain "Material", and fields passed
|
||||
/// will do 1 of 3 things:
|
||||
/// 1) They will not pass if they are an array associated with the whole mesh,
|
||||
/// 2) They will simply be passed to new cells if the array is associated with the cell set
|
||||
/// 3) They will be interpolated to new point locations if the array is associated with the point set
|
||||
///
|
||||
/// This algorithm requires passing a cell set of volume fraction information, not a point cell set.
|
||||
/// The exact fields are required:
|
||||
/// 1) A length cell set that specifies the number of materials associated to the cell.
|
||||
/// 2) A position cell set (or offset cell set) that specifies where the material IDs and VFs occur in the ID and VF arrays.
|
||||
/// 3) An ID array (whole array set) that stores the material ID information
|
||||
/// 4) An VF array (whole array set) that stores the fractional volume information for the respective material ID.
|
||||
/// Note that the cell VF information should add up to 1.0 across all materials for the cell, however this isn't checked in the code and might
|
||||
/// lead to undesirable results when iterating.
|
||||
///
|
||||
/// Note that this algorithm does not guarantee that the newly constructed cells will match the provided
|
||||
/// volume fractions, nor does it guarantee that there will exist a subcell of every material ID from the original cell.
|
||||
/// This usually occurs when the resolution of the mesh is too low (isolated materials in a single cell).
|
||||
///
|
||||
/// If wanted, this algorithm can iterate, adjusting cell VFs based on distance from the target values and the previous calculated iteration.
|
||||
/// This is done by setting the max iterations >0. In addition, the max percent error will allow for the filter to return early if the
|
||||
/// total error % of the entire dataset is less than the specified amount (defaults to 1.0, returns after first iteration). Finally,
|
||||
/// the error scaling and scaling decay allows for setting how much the cell VFs should react to the delta between target and calculated cell VFs.
|
||||
/// the error scaling will decay by the decay variable every iteration (multiplicitively).
|
||||
class MIRFilter : public vtkm::filter::FilterDataSet<MIRFilter>
|
||||
{
|
||||
public:
|
||||
/// @brief Sets the name of the offset/position cellset field in the dataset passed to the filter
|
||||
VTKM_CONT void SetPositionCellSetName(std::string name) { this->pos_name = name; }
|
||||
/// @brief Sets the name of the length cellset field in the dataset passed to the filter
|
||||
VTKM_CONT void SetLengthCellSetName(std::string name) { this->len_name = name; }
|
||||
/// @brief Sets the name of the ID whole-array set field in the dataset passed to the filter
|
||||
VTKM_CONT void SetIDWholeSetName(std::string name) { this->id_name = name; }
|
||||
/// @brief Sets the name of the VF whole-array set field in the dataset passed to the filter
|
||||
VTKM_CONT void SetVFWholeSetName(std::string name) { this->vf_name = name; }
|
||||
VTKM_CONT void SetMaxPercentError(vtkm::Float64 ma) { this->max_error = ma; }
|
||||
VTKM_CONT void SetMaxIterations(vtkm::IdComponent ma) { this->max_iter = ma; }
|
||||
VTKM_CONT void SetErrorScaling(vtkm::Float64 sc) { this->error_scaling = sc; }
|
||||
VTKM_CONT void SetScalingDecay(vtkm::Float64 sc) { this->scaling_decay = sc; }
|
||||
/// @brief Gets the output cell-set field name for the filter
|
||||
VTKM_CONT std::string GetOutputFieldName() { return this->OutputFieldName; }
|
||||
/// @brief Sets the output cell-set field name for the filter
|
||||
VTKM_CONT void SetOutputFieldName(std::string name) { this->OutputFieldName = name; }
|
||||
template <typename DerivedPolicy>
|
||||
VTKM_CONT vtkm::cont::DataSet DoExecute(const vtkm::cont::DataSet& input,
|
||||
vtkm::filter::PolicyBase<DerivedPolicy>);
|
||||
|
||||
template <typename T, typename StorageType, typename DerivedPolicy>
|
||||
VTKM_CONT bool DoMapField(vtkm::cont::DataSet& result,
|
||||
const vtkm::cont::ArrayHandle<T, StorageType>& input1,
|
||||
const vtkm::filter::FieldMetadata& fieldMeta,
|
||||
vtkm::filter::PolicyBase<DerivedPolicy>)
|
||||
{
|
||||
(void)result;
|
||||
(void)input1;
|
||||
(void)fieldMeta;
|
||||
|
||||
if (fieldMeta.GetName().compare(this->pos_name) == 0 ||
|
||||
fieldMeta.GetName().compare(this->len_name) == 0 ||
|
||||
fieldMeta.GetName().compare(this->id_name) == 0 ||
|
||||
fieldMeta.GetName().compare(this->vf_name) == 0)
|
||||
{
|
||||
// Remember, we will map the field manually...
|
||||
// Technically, this will be for all of them...thus ignore it
|
||||
return false;
|
||||
}
|
||||
vtkm::cont::ArrayHandle<T, vtkm::cont::StorageTagBasic> output;
|
||||
if (fieldMeta.IsPointField())
|
||||
{
|
||||
this->ProcessPointField1(input1, output);
|
||||
}
|
||||
else if (fieldMeta.IsCellField())
|
||||
{
|
||||
this->ProcessCellField(input1, output);
|
||||
}
|
||||
else
|
||||
{
|
||||
return false;
|
||||
}
|
||||
result.AddField(fieldMeta.AsField(output));
|
||||
return true;
|
||||
}
|
||||
|
||||
private:
|
||||
// Linear storage requirement, scales with size of point output
|
||||
template <typename T, typename StorageType, typename StorageType2>
|
||||
VTKM_CONT void ProcessPointField1(const vtkm::cont::ArrayHandle<T, StorageType>& input,
|
||||
vtkm::cont::ArrayHandle<T, StorageType2>& output);
|
||||
|
||||
// NOTE: The below assumes that MIR will not change the cell values when subdividing.
|
||||
template <typename T, typename StorageType, typename StorageType2>
|
||||
VTKM_CONT void ProcessCellField(const vtkm::cont::ArrayHandle<T, StorageType>& input,
|
||||
vtkm::cont::ArrayHandle<T, StorageType2>& output)
|
||||
{
|
||||
// Use a temporary permutation array to simplify the mapping:
|
||||
auto tmp = vtkm::cont::make_ArrayHandlePermutation(this->filterCellInterp, input);
|
||||
vtkm::cont::ArrayCopy(tmp, output);
|
||||
}
|
||||
|
||||
std::string pos_name;
|
||||
std::string len_name;
|
||||
std::string id_name;
|
||||
std::string vf_name;
|
||||
std::string OutputFieldName = std::string("cellMat");
|
||||
vtkm::Float64 max_error = vtkm::Float64(1.0);
|
||||
vtkm::Float64 scaling_decay = vtkm::Float64(1.0);
|
||||
vtkm::IdComponent max_iter = vtkm::IdComponent(0);
|
||||
vtkm::Float64 error_scaling = vtkm::Float64(0.0);
|
||||
vtkm::cont::ArrayHandle<vtkm::Id> filterCellInterp;
|
||||
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Float64, 8>> MIRWeights;
|
||||
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 8>> MIRIDs;
|
||||
};
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
#ifndef vtk_m_filter_MIRFilter_hxx
|
||||
#include <vtkm/filter/MIRFilter.hxx>
|
||||
#endif
|
||||
|
||||
#endif
|
||||
343
vtkm/filter/MIRFilter.hxx
Normal file
343
vtkm/filter/MIRFilter.hxx
Normal file
@ -0,0 +1,343 @@
|
||||
//============================================================================
|
||||
// Copyright (c) Kitware, Inc.
|
||||
// All rights reserved.
|
||||
// See LICENSE.txt for details.
|
||||
//
|
||||
// This software is distributed WITHOUT ANY WARRANTY; without even
|
||||
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
|
||||
// PURPOSE. See the above copyright notice for more information.
|
||||
//============================================================================
|
||||
|
||||
#ifndef vtk_m_filter_MIRFilter_hxx
|
||||
#define vtk_m_filter_MIRFilter_hxx
|
||||
|
||||
#include <vtkm/CellShape.h>
|
||||
#include <vtkm/Types.h>
|
||||
|
||||
#include <vtkm/cont/Algorithm.h>
|
||||
#include <vtkm/cont/ArrayCopy.h>
|
||||
#include <vtkm/cont/ArrayHandle.h>
|
||||
#include <vtkm/cont/ArrayHandlePermutation.h>
|
||||
#include <vtkm/cont/ArrayHandleView.h>
|
||||
#include <vtkm/cont/CellSetExplicit.h>
|
||||
#include <vtkm/cont/CellSetPermutation.h>
|
||||
#include <vtkm/cont/CoordinateSystem.h>
|
||||
#include <vtkm/cont/DynamicCellSet.h>
|
||||
#include <vtkm/cont/ErrorFilterExecution.h>
|
||||
#include <vtkm/cont/ExecutionObjectBase.h>
|
||||
#include <vtkm/cont/ImplicitFunctionHandle.h>
|
||||
#include <vtkm/cont/Timer.h>
|
||||
|
||||
#include <vtkm/exec/FunctorBase.h>
|
||||
|
||||
#include <vtkm/worklet/DispatcherMapField.h>
|
||||
#include <vtkm/worklet/DispatcherMapTopology.h>
|
||||
#include <vtkm/worklet/DispatcherReduceByKey.h>
|
||||
#include <vtkm/worklet/Keys.h>
|
||||
#include <vtkm/worklet/MIR.h>
|
||||
#include <vtkm/worklet/ScatterCounting.h>
|
||||
#include <vtkm/worklet/WorkletMapField.h>
|
||||
#include <vtkm/worklet/WorkletMapTopology.h>
|
||||
#include <vtkm/worklet/WorkletReduceByKey.h>
|
||||
#include <vtkm/worklet/clip/ClipTables.h>
|
||||
|
||||
#include <vtkm/filter/MeshQuality.h>
|
||||
|
||||
namespace vtkm
|
||||
{
|
||||
/*
|
||||
Todos:
|
||||
Enable some sort of cell culling, so if a cell doesn't have any work to do, it doesn't get called in future invocations of MIR
|
||||
|
||||
*/
|
||||
namespace filter
|
||||
{
|
||||
|
||||
template <typename T, typename StorageType, typename StorageType2>
|
||||
inline VTKM_CONT void MIRFilter::ProcessPointField1(
|
||||
const vtkm::cont::ArrayHandle<T, StorageType>& input,
|
||||
vtkm::cont::ArrayHandle<T, StorageType2>& output)
|
||||
{
|
||||
//VTKM_LOG_S(vtkm::cont::LogLevel::Warn, "HERE");
|
||||
vtkm::worklet::DestructPointWeightList destructWeightList;
|
||||
this->Invoke(destructWeightList, this->MIRIDs, this->MIRWeights, input, output);
|
||||
}
|
||||
//-----------------------------------------------------------------------------
|
||||
template <typename DerivedPolicy>
|
||||
inline VTKM_CONT vtkm::cont::DataSet MIRFilter::DoExecute(
|
||||
const vtkm::cont::DataSet& input,
|
||||
vtkm::filter::PolicyBase<DerivedPolicy> policy)
|
||||
{
|
||||
|
||||
|
||||
//{
|
||||
//(void)input;
|
||||
//(void)policy;
|
||||
vtkm::worklet::CheckFor2D cellCheck;
|
||||
vtkm::cont::ArrayHandle<vtkm::Id> count2D, count3D, countBad;
|
||||
this->Invoke(cellCheck,
|
||||
vtkm::filter::ApplyPolicyCellSet(input.GetCellSet(), policy, *this),
|
||||
count2D,
|
||||
count3D,
|
||||
countBad);
|
||||
vtkm::Id c2 = vtkm::cont::Algorithm::Reduce(count2D, vtkm::Id(0));
|
||||
vtkm::Id c3 = vtkm::cont::Algorithm::Reduce(count3D, vtkm::Id(0));
|
||||
vtkm::Id cB = vtkm::cont::Algorithm::Reduce(countBad, vtkm::Id(0));
|
||||
if (cB > vtkm::Id(0))
|
||||
{
|
||||
VTKM_LOG_S(
|
||||
vtkm::cont::LogLevel::Fatal,
|
||||
"Bad cell found in MIR filter input! Strictly only 2D -xor- 3D cell sets are permitted!");
|
||||
}
|
||||
if (c2 > vtkm::Id(0) && c3 > vtkm::Id(0))
|
||||
{
|
||||
VTKM_LOG_S(
|
||||
vtkm::cont::LogLevel::Fatal,
|
||||
"Bad cell mix found in MIR filter input! Input is not allowed to have both 2D and 3D cells.");
|
||||
}
|
||||
if (c2 == vtkm::Id(0) && c3 == vtkm::Id(0))
|
||||
{
|
||||
VTKM_LOG_S(vtkm::cont::LogLevel::Fatal,
|
||||
"No cells found for MIR filter! Please don't call me with nothing!");
|
||||
}
|
||||
|
||||
const vtkm::cont::CoordinateSystem inputCoords =
|
||||
input.GetCoordinateSystem(this->GetActiveCoordinateSystemIndex());
|
||||
vtkm::cont::ArrayHandle<vtkm::Float64> avgSizeTot;
|
||||
vtkm::worklet::MeshQuality<vtkm::filter::CellMetric> getVol;
|
||||
getVol.SetMetric(c3 > 0 ? vtkm::filter::CellMetric::VOLUME : vtkm::filter::CellMetric::AREA);
|
||||
this->Invoke(getVol,
|
||||
vtkm::filter::ApplyPolicyCellSet(input.GetCellSet(), policy, *this),
|
||||
inputCoords.GetData(),
|
||||
avgSizeTot);
|
||||
// First, load up all fields...
|
||||
vtkm::cont::Field or_pos = input.GetField(this->pos_name);
|
||||
vtkm::cont::Field or_len = input.GetField(this->len_name);
|
||||
vtkm::cont::Field or_ids = input.GetField(this->id_name);
|
||||
vtkm::cont::Field or_vfs = input.GetField(this->vf_name);
|
||||
// TODO: Check all fields for 'IsFieldCell'
|
||||
vtkm::cont::ArrayHandle<vtkm::Float32> vfsdata_or, vfsdata;
|
||||
vtkm::cont::ArrayHandle<vtkm::Id> idsdata_or, idsdata, lendata_or, lendata, posdata_or, posdata,
|
||||
allids;
|
||||
or_pos.GetData().AsArrayHandle(posdata_or);
|
||||
or_len.GetData().AsArrayHandle(lendata_or);
|
||||
or_ids.GetData().AsArrayHandle(idsdata_or);
|
||||
or_vfs.GetData().AsArrayHandle(vfsdata_or);
|
||||
vtkm::cont::ArrayCopy(idsdata_or, allids);
|
||||
vtkm::cont::Algorithm::Sort(allids);
|
||||
vtkm::cont::Algorithm::Unique(allids);
|
||||
vtkm::IdComponent numIDs = static_cast<vtkm::IdComponent>(allids.GetNumberOfValues());
|
||||
//using PortalConstType = vtkm::cont::ArrayHandle<vtkm::Id>::PortalConstControl;
|
||||
//PortalConstType readPortal = allids.GetPortalConstControl();
|
||||
using PortalConstType = vtkm::cont::ArrayHandle<vtkm::Id>::ReadPortalType;
|
||||
PortalConstType readPortal = allids.ReadPortal();
|
||||
vtkm::cont::ArrayCopy(idsdata_or, idsdata);
|
||||
vtkm::cont::ArrayCopy(lendata_or, lendata);
|
||||
vtkm::cont::ArrayCopy(posdata_or, posdata);
|
||||
vtkm::cont::ArrayCopy(vfsdata_or, vfsdata);
|
||||
//}
|
||||
|
||||
vtkm::cont::DataSet saved;
|
||||
// % error of the whole system, multiplied by the number of cells
|
||||
vtkm::Float64 totalError = this->max_error + vtkm::Float64(1.1); // Dummy value
|
||||
vtkm::IdComponent currentIterationNum = 0;
|
||||
|
||||
vtkm::worklet::MIRCases::MIRTables faceTableArray;
|
||||
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Float64, 8>> pointWeights;
|
||||
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 8>> pointIDs;
|
||||
vtkm::worklet::ConstructCellWeightList constructReverseInformation;
|
||||
vtkm::cont::ArrayHandleIndex pointCounter(input.GetNumberOfPoints());
|
||||
this->Invoke(constructReverseInformation, pointCounter, pointIDs, pointWeights);
|
||||
do
|
||||
{
|
||||
saved = vtkm::cont::DataSet();
|
||||
//auto cs = vtkm::filter::ApplyPolicyCellSet(input.GetCellSet(), policy, *this);
|
||||
//saved.SetCellSet(cs);
|
||||
saved.SetCellSet(input.GetCellSet());
|
||||
vtkm::cont::CoordinateSystem outCo(inputCoords.GetName(), inputCoords.GetData());
|
||||
saved.AddCoordinateSystem(outCo);
|
||||
//{
|
||||
vtkm::cont::ArrayHandle<vtkm::Id> currentcellIDs;
|
||||
vtkm::cont::ArrayHandle<vtkm::Id> pointlen, pointpos, pointid;
|
||||
vtkm::cont::ArrayHandle<vtkm::Float64> pointvf;
|
||||
vtkm::worklet::CombineVFsForPoints_C convertOrigCellTo;
|
||||
vtkm::worklet::CombineVFsForPoints convertOrigCellTo_Full;
|
||||
|
||||
|
||||
this->Invoke(convertOrigCellTo, saved.GetCellSet(), lendata, posdata, idsdata, pointlen);
|
||||
vtkm::Id pointcount = vtkm::cont::Algorithm::ScanExclusive(pointlen, pointpos);
|
||||
pointvf.Allocate(pointcount);
|
||||
pointid.Allocate(pointcount);
|
||||
this->Invoke(convertOrigCellTo_Full,
|
||||
saved.GetCellSet(),
|
||||
lendata,
|
||||
posdata,
|
||||
idsdata,
|
||||
vfsdata,
|
||||
pointpos,
|
||||
pointid,
|
||||
pointvf);
|
||||
|
||||
|
||||
//}
|
||||
vtkm::worklet::MIRObject<vtkm::Id, vtkm::Float64> mirobj(
|
||||
pointlen, pointpos, pointid, pointvf); // This is point VF data...
|
||||
auto tmp2 = std::vector<vtkm::Id>(saved.GetCellSet().GetNumberOfCells(), -1);
|
||||
vtkm::cont::ArrayHandle<vtkm::Id> prevMat =
|
||||
vtkm::cont::make_ArrayHandle(tmp2, vtkm::CopyFlag::On);
|
||||
vtkm::cont::ArrayHandle<vtkm::Id> cellLookback;
|
||||
vtkm::cont::ArrayHandleIndex tmp_ind(saved.GetCellSet().GetNumberOfCells());
|
||||
vtkm::cont::ArrayCopy(tmp_ind, cellLookback);
|
||||
vtkm::IdComponent currentMatLoc = 0;
|
||||
|
||||
|
||||
|
||||
|
||||
while (currentMatLoc < numIDs)
|
||||
{
|
||||
vtkm::IdComponent currentMatID =
|
||||
static_cast<vtkm::IdComponent>(readPortal.Get(currentMatLoc++));
|
||||
if (currentMatID < 1)
|
||||
{
|
||||
VTKM_LOG_S(
|
||||
vtkm::cont::LogLevel::Fatal,
|
||||
"MIR filter does not accept materials with an non-positive ID! Material id in offense: "
|
||||
<< currentMatID
|
||||
<< ". Please remap all ID values to only positive numbers to avoid this issue.");
|
||||
}
|
||||
// First go through and pick out the previous and current material VFs for each cell.
|
||||
//{
|
||||
vtkm::worklet::ExtractVFsForMIR_C extractCurrentMatVF;
|
||||
vtkm::cont::ArrayHandle<vtkm::Id> currentCellPointCounts;
|
||||
this->Invoke(extractCurrentMatVF, saved.GetCellSet(), currentCellPointCounts);
|
||||
vtkm::worklet::ExtractVFsForMIR extractCurrentMatVF_SC(currentMatID);
|
||||
vtkm::worklet::ScatterCounting extractCurrentMatVF_SC_scatter =
|
||||
extractCurrentMatVF_SC.MakeScatter(currentCellPointCounts);
|
||||
vtkm::cont::ArrayHandle<vtkm::Float64> currentMatVF;
|
||||
vtkm::cont::ArrayHandle<vtkm::Float64> previousMatVF;
|
||||
this->Invoke(extractCurrentMatVF_SC,
|
||||
extractCurrentMatVF_SC_scatter,
|
||||
saved.GetCellSet(),
|
||||
mirobj,
|
||||
prevMat,
|
||||
currentMatVF,
|
||||
previousMatVF);
|
||||
//}
|
||||
// Next see if we need to perform any work at all...
|
||||
if (currentMatLoc != 0)
|
||||
{
|
||||
// Run MIR, possibly changing colors...
|
||||
vtkm::cont::ArrayHandle<vtkm::Id> cellVFPointOffsets;
|
||||
vtkm::cont::Algorithm::ScanExclusive(currentCellPointCounts, cellVFPointOffsets);
|
||||
vtkm::worklet::MIR mir;
|
||||
vtkm::cont::ArrayHandle<vtkm::Id> newCellLookback, newCellID;
|
||||
|
||||
|
||||
vtkm::cont::CellSetExplicit<> out = mir.Run(saved.GetCellSet(),
|
||||
previousMatVF,
|
||||
currentMatVF,
|
||||
cellVFPointOffsets,
|
||||
prevMat,
|
||||
currentMatID,
|
||||
cellLookback,
|
||||
newCellID,
|
||||
newCellLookback);
|
||||
vtkm::cont::ArrayCopy(newCellLookback, cellLookback);
|
||||
vtkm::cont::ArrayCopy(newCellID, prevMat);
|
||||
auto data = saved.GetCoordinateSystem(0).GetDataAsMultiplexer();
|
||||
auto coords = mir.ProcessPointField(data);
|
||||
// Now convert the point VFs...
|
||||
vtkm::cont::ArrayHandle<vtkm::Id> plen, ppos, pids;
|
||||
vtkm::cont::ArrayHandle<vtkm::Float64> pvf;
|
||||
mir.ProcessMIRField(mirobj.getPointLenArr(),
|
||||
mirobj.getPointPosArr(),
|
||||
mirobj.getPointIDArr(),
|
||||
mirobj.getPointVFArr(),
|
||||
plen,
|
||||
ppos,
|
||||
pids,
|
||||
pvf);
|
||||
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Float64, 8>> tmppointWeights;
|
||||
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 8>> tmppointIDs;
|
||||
mir.ProcessSimpleMIRField(pointIDs, pointWeights, tmppointIDs, tmppointWeights);
|
||||
vtkm::cont::ArrayCopy(tmppointIDs, pointIDs);
|
||||
vtkm::cont::ArrayCopy(tmppointWeights, pointWeights);
|
||||
//FileSaver fs;
|
||||
//fs(("pID" + std::to_string(currentMatID) + ".txt").c_str(), pointIDs);
|
||||
//fs(("wID" + std::to_string(currentMatID) + ".txt").c_str(), pointWeights);
|
||||
mirobj = vtkm::worklet::MIRObject<vtkm::Id, vtkm::Float64>(plen, ppos, pids, pvf);
|
||||
saved = vtkm::cont::DataSet();
|
||||
saved.SetCellSet(out);
|
||||
vtkm::cont::CoordinateSystem outCo2(inputCoords.GetName(), coords);
|
||||
saved.AddCoordinateSystem(outCo2);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// Hacking workaround to not clone an entire dataset.
|
||||
vtkm::cont::ArrayHandle<vtkm::Float64> avgSize;
|
||||
this->Invoke(getVol, saved.GetCellSet(), saved.GetCoordinateSystem(0).GetData(), avgSize);
|
||||
|
||||
vtkm::worklet::CalcError_C calcErrC;
|
||||
vtkm::worklet::Keys<vtkm::Id> cellKeys(cellLookback);
|
||||
vtkm::cont::ArrayCopy(cellLookback, filterCellInterp);
|
||||
vtkm::cont::ArrayHandle<vtkm::Id> lenOut, posOut, idsOut;
|
||||
vtkm::cont::ArrayHandle<vtkm::Float64> vfsOut, totalErrorOut;
|
||||
|
||||
lenOut.Allocate(cellKeys.GetUniqueKeys().GetNumberOfValues());
|
||||
this->Invoke(calcErrC, cellKeys, prevMat, lendata_or, posdata_or, idsdata_or, lenOut);
|
||||
|
||||
vtkm::Id numIDsOut = vtkm::cont::Algorithm::ScanExclusive(lenOut, posOut);
|
||||
idsOut.Allocate(numIDsOut);
|
||||
vfsOut.Allocate(numIDsOut);
|
||||
vtkm::worklet::CalcError calcErr(this->error_scaling);
|
||||
this->Invoke(calcErr,
|
||||
cellKeys,
|
||||
prevMat,
|
||||
avgSize,
|
||||
lendata_or,
|
||||
posdata_or,
|
||||
idsdata_or,
|
||||
vfsdata_or,
|
||||
lendata,
|
||||
posdata,
|
||||
idsdata,
|
||||
vfsdata,
|
||||
lenOut,
|
||||
posOut,
|
||||
idsOut,
|
||||
vfsOut,
|
||||
avgSizeTot,
|
||||
totalErrorOut);
|
||||
totalError = vtkm::cont::Algorithm::Reduce(totalErrorOut, vtkm::Float64(0));
|
||||
vtkm::cont::ArrayCopy(lenOut, lendata);
|
||||
vtkm::cont::ArrayCopy(posOut, posdata);
|
||||
vtkm::cont::ArrayCopy(idsOut, idsdata);
|
||||
vtkm::cont::ArrayCopy(vfsOut, vfsdata);
|
||||
// Clean up the cells by calculating their volumes, and then calculate the relative error for each cell.
|
||||
// Note that the total error needs to be rescaled by the number of cells to get the % error.
|
||||
totalError =
|
||||
totalError /
|
||||
vtkm::Float64(
|
||||
vtkm::filter::ApplyPolicyCellSet(input.GetCellSet(), policy, *this).GetNumberOfCells());
|
||||
this->error_scaling *= this->scaling_decay;
|
||||
|
||||
|
||||
VTKM_LOG_S(vtkm::cont::LogLevel::Info,
|
||||
"Mir iteration " << currentIterationNum + 1 << "/" << this->max_iter
|
||||
<< "\t Total error: " << totalError);
|
||||
saved.AddField(vtkm::cont::Field(
|
||||
this->GetOutputFieldName(), vtkm::cont::Field::Association::CELL_SET, prevMat));
|
||||
|
||||
vtkm::cont::ArrayCopy(pointIDs, this->MIRIDs);
|
||||
vtkm::cont::ArrayCopy(pointWeights, this->MIRWeights);
|
||||
} while ((++currentIterationNum <= this->max_iter) && totalError >= this->max_error);
|
||||
|
||||
|
||||
return saved;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
#endif
|
||||
@ -50,6 +50,7 @@ set(unit_tests
|
||||
UnitTestMaskFilter.cxx
|
||||
UnitTestMaskPointsFilter.cxx
|
||||
UnitTestMeshQualityFilter.cxx
|
||||
UnitTestMIRFilter.cxx
|
||||
UnitTestMultiBlockFilter.cxx
|
||||
UnitTestNDEntropyFilter.cxx
|
||||
UnitTestNDHistogramFilter.cxx
|
||||
|
||||
123
vtkm/filter/testing/UnitTestMIRFilter.cxx
Normal file
123
vtkm/filter/testing/UnitTestMIRFilter.cxx
Normal file
@ -0,0 +1,123 @@
|
||||
//============================================================================
|
||||
// Copyright (c) Kitware, Inc.
|
||||
// All rights reserved.
|
||||
// See LICENSE.txt for details.
|
||||
//
|
||||
// This software is distributed WITHOUT ANY WARRANTY; without even
|
||||
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
|
||||
// PURPOSE. See the above copyright notice for more information.
|
||||
//============================================================================
|
||||
|
||||
#include <vtkm/cont/DataSetBuilderExplicit.h>
|
||||
#include <vtkm/cont/testing/Testing.h>
|
||||
|
||||
#include <vtkm/filter/MIRFilter.h>
|
||||
|
||||
void ConnectionHelperHex(std::vector<vtkm::Id>& conn, int x, int y, int z, int mx, int my, int mz)
|
||||
{
|
||||
(void)mz;
|
||||
conn.push_back(mx * (my * z + y) + x);
|
||||
conn.push_back(mx * (my * z + y) + x + 1);
|
||||
conn.push_back(mx * (my * z + y + 1) + x + 1);
|
||||
conn.push_back(mx * (my * z + y + 1) + x);
|
||||
conn.push_back(mx * (my * (z + 1) + y) + x);
|
||||
conn.push_back(mx * (my * (z + 1) + y) + x + 1);
|
||||
conn.push_back(mx * (my * (z + 1) + y + 1) + x + 1);
|
||||
conn.push_back(mx * (my * (z + 1) + y + 1) + x);
|
||||
}
|
||||
|
||||
vtkm::cont::DataSet GetTestDataSet()
|
||||
{
|
||||
vtkm::cont::DataSetBuilderExplicit dsb;
|
||||
|
||||
int mx = 3, my = 3, mz = 3;
|
||||
|
||||
|
||||
std::vector<vtkm::UInt8> shapes;
|
||||
std::vector<vtkm::Id> connections;
|
||||
std::vector<vtkm::IdComponent> numberofInd;
|
||||
std::vector<vtkm::Vec3f_32> points;
|
||||
|
||||
for (int z = 0; z < mz - 1; z++)
|
||||
{
|
||||
for (int y = 0; y < my - 1; y++)
|
||||
{
|
||||
for (int x = 0; x < mx - 1; x++)
|
||||
{
|
||||
ConnectionHelperHex(connections, x, y, z, mx, my, mz);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
std::vector<vtkm::Id> idAR{ 1, 2, 2, 1, 2, 1, 1, 2 };
|
||||
std::vector<vtkm::Id> lnAR{ 1, 1, 1, 1, 1, 1, 1, 1 };
|
||||
std::vector<vtkm::Id> ofAR{ 0, 1, 2, 3, 4, 5, 6, 7 };
|
||||
vtkm::cont::ArrayHandle<vtkm::Id> offsets =
|
||||
vtkm::cont::make_ArrayHandle(ofAR, vtkm::CopyFlag::On);
|
||||
vtkm::cont::ArrayHandle<vtkm::Id> lengths =
|
||||
vtkm::cont::make_ArrayHandle(lnAR, vtkm::CopyFlag::On);
|
||||
vtkm::cont::ArrayHandle<vtkm::Id> ids = vtkm::cont::make_ArrayHandle(idAR, vtkm::CopyFlag::On);
|
||||
std::vector<vtkm::Float32> vfAR{ 1, 1, 1, 1, 1, 1, 1, 1 };
|
||||
vtkm::cont::ArrayHandle<vtkm::Float32> vfs =
|
||||
vtkm::cont::make_ArrayHandle(vfAR, vtkm::CopyFlag::On);
|
||||
|
||||
shapes.reserve((mx - 1) * (my - 1) * (mz - 1));
|
||||
numberofInd.reserve((mx - 1) * (my - 1) * (mz - 1));
|
||||
for (int i = 0; i < (mx - 1) * (my - 1) * (mz - 1); i++)
|
||||
{
|
||||
shapes.push_back(vtkm::CELL_SHAPE_HEXAHEDRON);
|
||||
numberofInd.push_back(8);
|
||||
}
|
||||
|
||||
points.reserve(mz * my * mx);
|
||||
for (int z = 0; z < mz; z++)
|
||||
{
|
||||
for (int y = 0; y < my; y++)
|
||||
{
|
||||
for (int x = 0; x < mx; x++)
|
||||
{
|
||||
vtkm::Vec3f_32 point(static_cast<vtkm::Float32>(x),
|
||||
static_cast<vtkm::Float32>(y),
|
||||
static_cast<vtkm::Float32>(z));
|
||||
points.push_back(point);
|
||||
}
|
||||
}
|
||||
}
|
||||
VTKM_LOG_S(vtkm::cont::LogLevel::Warn, "A");
|
||||
vtkm::cont::DataSet ds = dsb.Create(points, shapes, numberofInd, connections);
|
||||
VTKM_LOG_S(vtkm::cont::LogLevel::Warn, "B");
|
||||
ds.AddField(vtkm::cont::Field("scatter_pos", vtkm::cont::Field::Association::CELL_SET, offsets));
|
||||
VTKM_LOG_S(vtkm::cont::LogLevel::Warn, "C");
|
||||
ds.AddField(vtkm::cont::Field("scatter_len", vtkm::cont::Field::Association::CELL_SET, lengths));
|
||||
ds.AddField(vtkm::cont::Field("scatter_ids", vtkm::cont::Field::Association::WHOLE_MESH, ids));
|
||||
ds.AddField(vtkm::cont::Field("scatter_vfs", vtkm::cont::Field::Association::WHOLE_MESH, vfs));
|
||||
|
||||
return ds;
|
||||
}
|
||||
|
||||
void TestMIR()
|
||||
{
|
||||
vtkm::cont::DataSet ds = GetTestDataSet();
|
||||
|
||||
vtkm::filter::MIRFilter mir;
|
||||
mir.SetIDWholeSetName("scatter_ids");
|
||||
mir.SetPositionCellSetName("scatter_pos");
|
||||
mir.SetLengthCellSetName("scatter_len");
|
||||
mir.SetVFWholeSetName("scatter_vfs");
|
||||
|
||||
mir.SetErrorScaling(vtkm::Float64(0.2));
|
||||
mir.SetScalingDecay(vtkm::Float64(1.0));
|
||||
mir.SetMaxIterations(vtkm::IdComponent(0)); // =0 -> No iterations..
|
||||
mir.SetMaxPercentError(vtkm::Float64(
|
||||
0.00001)); // Only useful for iterations >= 1, will stop iterating if total % error for entire mesh is less than this value
|
||||
// Note it is mathematically impossible to obtain 0% error outside of VERY special cases (neglecting float error)
|
||||
vtkm::cont::DataSet ds_from_mir = mir.Execute(ds);
|
||||
|
||||
// Test if ds_from_mir has 40 cells
|
||||
VTKM_TEST_ASSERT(ds_from_mir.GetNumberOfCells() == 40, "Wrong number of output cells");
|
||||
}
|
||||
|
||||
int UnitTestMIRFilter(int argc, char* argv[])
|
||||
{
|
||||
return vtkm::cont::testing::Testing::Run(TestMIR, argc, argv);
|
||||
}
|
||||
@ -48,6 +48,7 @@ set(headers
|
||||
MaskPoints.h
|
||||
MaskSelect.h
|
||||
MeshQuality.h
|
||||
MIR.h
|
||||
NDimsEntropy.h
|
||||
NDimsHistMarginalization.h
|
||||
NDimsHistogram.h
|
||||
@ -136,6 +137,7 @@ add_subdirectory(cosmotools)
|
||||
add_subdirectory(gradient)
|
||||
add_subdirectory(histogram)
|
||||
add_subdirectory(lcs)
|
||||
add_subdirectory(mir)
|
||||
add_subdirectory(moments)
|
||||
add_subdirectory(splatkernels)
|
||||
add_subdirectory(spatialstructure)
|
||||
|
||||
2428
vtkm/worklet/MIR.h
Normal file
2428
vtkm/worklet/MIR.h
Normal file
File diff suppressed because it is too large
Load Diff
15
vtkm/worklet/mir/CMakeLists.txt
Normal file
15
vtkm/worklet/mir/CMakeLists.txt
Normal file
@ -0,0 +1,15 @@
|
||||
##============================================================================
|
||||
## Copyright (c) Kitware, Inc.
|
||||
## All rights reserved.
|
||||
## See LICENSE.txt for details.
|
||||
##
|
||||
## This software is distributed WITHOUT ANY WARRANTY; without even
|
||||
## the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
|
||||
## PURPOSE. See the above copyright notice for more information.
|
||||
##============================================================================
|
||||
|
||||
set(headers
|
||||
MIRTables.h
|
||||
)
|
||||
|
||||
vtkm_declare_headers(${headers})
|
||||
11457
vtkm/worklet/mir/MIRTables.h
Normal file
11457
vtkm/worklet/mir/MIRTables.h
Normal file
File diff suppressed because it is too large
Load Diff
Loading…
Reference in New Issue
Block a user