Commit 9989fad0 authored by js9's avatar js9

Adding material iteration capability, primarily for Zoo algorithm.



git-svn-id: http://visit.ilight.com/svn/visit/trunk/src@6285 18c085ea-50e0-402c-830e-de6fd14e8384
parent ff6325ab
......@@ -126,7 +126,7 @@ class DBATTS_API avtSILMatrix
SILCategoryRole GetRole1(void) const { return role1; };
SILCategoryRole GetRole2(void) const { return role2; };
SILCategoryRole GetRoleForCollection(int ind) const
{ return (ind < set1.size() ? role2 : role1); };
{ return (ind < int(set1.size()) ? role2 : role1); };
protected:
avtSIL *sil;
......
......@@ -3461,6 +3461,9 @@ avtGenericDatabase::AddOriginalNodesArray(vtkDataSet *ds, const int domain)
// they had only one material. Added code to account for that in
// case of creating boundary surfaces.
//
// Jeremy Meredith, Fri Feb 13 11:22:39 EST 2009
// Added MIR iteration capability.
//
// ****************************************************************************
avtDataTree_p
......@@ -3476,6 +3479,8 @@ avtGenericDatabase::MaterialSelect(vtkDataSet *ds, avtMaterial *mat,
bool simplifyHeavilyMixedZones,
int maxMatsPerZone,
int mirAlgorithm,
int mirNumIterations,
float mirIterationDamping,
float isovolumeMIRVF,
bool didGhosts,
bool &subdivisionOccurred,
......@@ -3535,7 +3540,9 @@ avtGenericDatabase::MaterialSelect(vtkDataSet *ds, avtMaterial *mat,
needValidConnectivity,
needSmoothMaterialInterfaces,
needCleanZonesOnly, simplifyHeavilyMixedZones,
maxMatsPerZone, mirAlgorithm, isovolumeMIRVF,
maxMatsPerZone, mirAlgorithm,
mirNumIterations, mirIterationDamping,
isovolumeMIRVF,
didGhosts,
subdivisionOccurred,
notAllCellsSubdivided, reUseMIR,
......@@ -4328,6 +4335,9 @@ avtGenericDatabase::SpeciesSelect(avtDatasetCollection &dsc,
// Hank Childs, Thu Feb 21 16:41:25 PST 2008
// Throw an exception if we get a bad MIR type value.
//
// Jeremy Meredith, Fri Feb 13 11:22:39 EST 2009
// Added MIR iteration capability.
//
// ****************************************************************************
void_ref_ptr
......@@ -4338,6 +4348,8 @@ avtGenericDatabase::GetMIR(int domain, const char *varname, int timestep,
bool needCleanZonesOnly,
bool simplifyHeavilyMixedZones, int maxMatsPerZone,
int mirAlgorithm,
int mirNumIterations,
float mirIterationDamping,
float isovolumeMIRVF,
bool didGhosts,
bool &subdivisionOccurred,
......@@ -4347,7 +4359,7 @@ avtGenericDatabase::GetMIR(int domain, const char *varname, int timestep,
mat_to_use = mat;
char cacheLbl[1000];
sprintf(cacheLbl, "MIR_%s_%s_%s_%s_%s_%d_%f_%s",
sprintf(cacheLbl, "MIR_%s_%s_%s_%s_%s_%d_%f_%s_%d_%f",
needValidConnectivity ? "FullSubdiv" : "MinimalSubdiv",
needSmoothMaterialInterfaces ? "Smooth" : "NotSmooth",
needCleanZonesOnly ? "CleanOnly" : "SplitMixed",
......@@ -4355,7 +4367,10 @@ avtGenericDatabase::GetMIR(int domain, const char *varname, int timestep,
simplifyHeavilyMixedZones ? "Simplify" : "NoSimplify",
maxMatsPerZone,
isovolumeMIRVF,
mirAlgorithm==0 ? "TetMIR" : (mirAlgorithm==1 ? "ZooMIR" : "IsovolumeMIR"));
mirAlgorithm==0 ? "TetMIR" : (mirAlgorithm==1 ? "ZooMIR" : "IsovolumeMIR"),
mirNumIterations,
mirIterationDamping);
//
// See if we already have the data lying around.
......@@ -4402,6 +4417,8 @@ avtGenericDatabase::GetMIR(int domain, const char *varname, int timestep,
}
mir->SetAlgorithm(mirAlgorithm);
mir->SetNumIterations(mirNumIterations);
mir->SetIterationDamping(mirIterationDamping);
mir->SetLeaveCleanZonesWhole(!needValidConnectivity);
mir->SetSmoothing(needSmoothMaterialInterfaces);
mir->SetCleanZonesOnly(needCleanZonesOnly);
......@@ -7195,6 +7212,9 @@ avtGenericDatabase::ApplyGhostForDomainNesting(avtDatasetCollection &ds,
// Hank Childs, Wed Jul 25 14:16:36 PDT 2007
// Renamed method: NeedBoundarySurfaces -> GetBoundarySurfaceRepresentation.
//
// Jeremy Meredith, Fri Feb 13 11:22:39 EST 2009
// Added MIR iteration capability.
//
// ****************************************************************************
void
......@@ -7269,6 +7289,8 @@ avtGenericDatabase::MaterialSelect(avtDatasetCollection &ds,
spec->SimplifyHeavilyMixedZones(),
spec->MaxMaterialsPerZone(),
spec->MIRAlgorithm(),
spec->MIRNumIterations(),
spec->MIRIterationDamping(),
spec->IsovolumeMIRVF(),
didGhosts, so, nacs, reUseMIR);
......
......@@ -330,6 +330,9 @@ class vtkUnstructuredGrid;
// are streaming, not about whether we are doing dynamic load balancing.
// And the two are no longer synonymous.
//
// Jeremy Meredith, Fri Feb 13 12:04:16 EST 2009
// Added MIR iteration capability.
//
// ****************************************************************************
class DATABASE_API avtGenericDatabase : public avtDatasetDatabase
......@@ -450,11 +453,13 @@ class DATABASE_API avtGenericDatabase : public avtDatasetDatabase
stringVector &,
stringVector &,
bool, bool, bool, bool, bool,
bool, bool, int, int, float, bool,
bool, bool, int, int,
int, float, float, bool,
bool&, bool&, bool);
void_ref_ptr GetMIR(int, const char *, int, vtkDataSet*,
avtMaterial *, int, bool, bool, bool,
bool, int, int, float, bool, bool&, bool&,bool,
bool, int, int, int, float,
float, bool, bool&, bool&,bool,
avtMaterial *&);
avtMaterial *GetMaterial(int, const char *, int, const avtDataRequest_p = 0);
avtSpecies *GetSpecies(int, const char *, int);
......
......@@ -119,6 +119,22 @@ MIR::SetNumIterations(int ni)
options.numIterations = ni;
}
// ****************************************************************************
// Method: MIR::SetIterationDamping
//
// Purpose:
// Set the option iterationDamping.
//
// Programmer: Jeremy Meredith
// Creation: February 12, 2009
//
// ****************************************************************************
void
MIR::SetIterationDamping(float damp)
{
options.iterationDamping = damp;
}
// ****************************************************************************
// Method: MIR::SetSmoothing
//
......
......@@ -78,6 +78,9 @@ class avtSpecies;
// Jeremy Meredith, Thu Aug 18 16:35:05 PDT 2005
// Added algorithm selector, and added VF for isovolume algorithm.
//
// Jeremy Meredith, Fri Feb 13 11:22:39 EST 2009
// Added MIR iteration capability.
//
// ****************************************************************************
class MIR_API MIR
{
......@@ -91,6 +94,7 @@ class MIR_API MIR
void SetAlgorithm(int);
void SetSubdivisionLevel(MIROptions::SubdivisionLevel);
void SetNumIterations(int);
void SetIterationDamping(float);
void SetSmoothing(bool);
void SetLeaveCleanZonesWhole(bool);
void SetCleanZonesOnly(bool);
......
......@@ -48,12 +48,16 @@
// Jeremy Meredith, Thu Aug 18 16:36:55 PDT 2005
// Added algorithm and isovolumeVF.
//
// Jeremy Meredith, Fri Feb 13 11:22:39 EST 2009
// Added MIR iteration capability.
//
// ****************************************************************************
MIROptions::MIROptions()
{
algorithm = 1;
subdivisionLevel = MIROptions::Low;
numIterations = 0;
iterationDamping = 0.4;
smoothing = true;
leaveCleanZonesWhole = true;
cleanZonesOnly = false;
......
......@@ -75,6 +75,9 @@
// Jeremy Meredith, Thu Aug 18 16:36:42 PDT 2005
// Added algorithm and isovolumeVF.
//
// Jeremy Meredith, Fri Feb 13 11:22:39 EST 2009
// Added MIR iteration and damping capability.
//
// ****************************************************************************
class MIR_API MIROptions
{
......@@ -89,6 +92,7 @@ class MIR_API MIROptions
int algorithm;
SubdivisionLevel subdivisionLevel;
int numIterations;
float iterationDamping;
bool smoothing;
bool leaveCleanZonesWhole;
bool cleanZonesOnly;
......
......@@ -97,14 +97,14 @@ class QuadraticHash
private:
unsigned int (*hashfunc)(K&);
unsigned int size;
int size;
int sizeindex;
unsigned int currentcell;
K currentKey;
int iteratorpos;
unsigned int numvalid;
int numvalid;
Entry *table;
};
......@@ -178,7 +178,7 @@ QuadraticHash<K,V>::~QuadraticHash()
template <class K, class V>
bool QuadraticHash<K,V>::Find(K &key)
{
unsigned int emptycell = -1;
int emptycell = -1;
unsigned int startcell = hashfunc(key) % size;;
unsigned int probeindex = 0;
......
......@@ -987,7 +987,7 @@ TetMIR::GetDataset(vector<int> mats, vtkDataSet *ds,
for (size_t i = 0; i < mats.size(); i++)
{
int origmatno = mats[i];
if (origmatno < mapMatToUsedMat.size())
if (origmatno < int(mapMatToUsedMat.size()))
{
int usedmatno = mapMatToUsedMat[origmatno];
if (usedmatno != -1)
......
......@@ -48,6 +48,8 @@
#include "BitUtils.h"
#include "ResampledMat.h"
#include <verdict.h>
// ****************************************************************************
// Constructor: CellReconstructor::CellReconstructor
//
......@@ -422,3 +424,110 @@ CellReconstructor::CreateOutputShape(TempCell &old,
outlist.push_back(cell);
}
// ****************************************************************************
// Method: CellReconstructor::CalculateVolumeOrAreaHelper
//
// Purpose:
// Calculate the area or volume of a VTK cell. Logic taken
// from the mesh quality expressions.
//
// Arguments:
// type vtk cell type
// coords the coordinates in VTK order
//
// Programmer: Jeremy Meredith
// Creation: February 13, 2009
//
// ****************************************************************************
inline void
Swap1(double &a, double &b)
{
double tmp = a;
a = b;
b = tmp;
}
inline void
Swap3(double c[][3], int a, int b)
{
Swap1(c[a][0], c[b][0]);
Swap1(c[a][1], c[b][1]);
Swap1(c[a][2], c[b][2]);
}
inline
void Copy3(double coords[][3], double a[], int i)
{
a[0] = coords[i][0];
a[1] = coords[i][1];
a[2] = coords[i][2];
}
double
CellReconstructor::CalculateVolumeOrAreaHelper(int type, double coords[][3])
{
switch (type)
{
case VTK_TRIANGLE:
return v_tri_area(3, coords);
case VTK_QUAD:
return v_quad_area(4, coords);
case VTK_PIXEL:
Swap3(coords, 2, 3);
return v_quad_area(4, coords);
case VTK_VOXEL:
Swap3(coords, 2,3);
Swap3(coords, 6,7);
return v_hex_volume(8,coords);
case VTK_HEXAHEDRON:
return v_hex_volume(8,coords);
case VTK_TETRA:
return v_tet_volume(4,coords);
case VTK_WEDGE:
{
int subdiv[3][4] = { {0,5,4,3}, {0,2,1,4}, {0,4,5,2} };
double tet_coords[4][3];
double vol = 0;
for (int i = 0 ; i < 3 ; i++)
{
for (int j = 0 ; j < 4 ; j++)
for (int k = 0 ; k < 3 ; k++)
tet_coords[j][k] = coords[subdiv[i][j]][k];
vol += v_tet_volume(4, tet_coords);
}
return vol;
}
// The verdict metric for pyramid I have yet to figure out how to work.
// However, it does the same thing that we do here: Divide the pyramid
// into two tetrahedrons.
case VTK_PYRAMID:
{
double one[4][3];
double two[4][3];
Copy3(coords,one[0], 0);
Copy3(coords,one[1], 1);
Copy3(coords,one[2], 2);
Copy3(coords,one[3], 4);
Copy3(coords,two[0], 0);
Copy3(coords,two[1], 2);
Copy3(coords,two[2], 3);
Copy3(coords,two[3], 4);
return v_tet_volume(4,one) + v_tet_volume(4,two);
}
}
// ERROR; SHOULDN'T GET HERE!!!!
return -99999;
}
......@@ -61,6 +61,12 @@
// Mark C. Miller, Thu Feb 9 21:06:10 PST 2006
// Renamed Array class to VisItArray to avoid name collisions with
// third-party libs
//
// Jeremy Meredith, Fri Feb 13 10:56:43 EST 2009
// Allowed ReconstructCell to take a new argument where it will, if
// desired, output the volume fractions for the reconstructed materials.
// Also, added helper function to calculate volume or area.
//
// ****************************************************************************
class CellReconstructor
{
......@@ -68,17 +74,21 @@ class CellReconstructor
CellReconstructor(vtkDataSet*, avtMaterial*, ResampledMat&, int, int, bool,
MIRConnectivity&, ZooMIR&);
virtual ~CellReconstructor();
virtual void ReconstructCell(int, int, int, int*) = 0;
virtual void ReconstructCell(int, int, int, int*, double*) = 0;
protected:
vtkDataSet *mesh;
avtMaterial *mat;
avtMaterial *origMat;
ResampledMat &rm;
int nPoints;
int nCells;
MIRConnectivity &conn;
ZooMIR &mir;
int nMaterials;
static double CalculateVolumeOrAreaHelper(int celltype,double coords[][3]);
protected:
int cellid;
int celltype;
......
......@@ -85,10 +85,14 @@ IsovolumeCellReconstructor::IsovolumeCellReconstructor(vtkDataSet *d,
// Jeremy Meredith, Thu Aug 18 17:09:36 PDT 2005
// Initial version. Lifted heavily from the Recursive version.
//
// Jeremy Meredith, Fri Feb 13 11:05:06 EST 2009
// Added calculation of output vf's per material, if requested.
//
// ****************************************************************************
void
IsovolumeCellReconstructor::ReconstructCell(int cellid_, int celltype_,
int nids_, int *ids_)
int nids_, int *ids_,
double *outputvfs)
{
cellid = cellid_;
celltype = celltype_;
......@@ -348,6 +352,17 @@ IsovolumeCellReconstructor::ReconstructCell(int cellid_, int celltype_,
}
}
// If we're going to calculate actual volume fractions, first
// zero them out, then accumulate the output cell partial contributions.
double totalvol = 0;
if (outputvfs)
{
for (int matno=0; matno < nMaterials; matno++)
{
outputvfs[matno] = 0.0;
}
}
// Spit the reconstructed cells into the output zone list
for (int out = 0 ; out < outlist.size() ; out++)
{
......@@ -364,6 +379,39 @@ IsovolumeCellReconstructor::ReconstructCell(int cellid_, int celltype_,
for (int n=0; n<outcell.nnodes; n++)
mir.indexList.push_back(outcell.ids[n]);
if (outputvfs)
{
double coords[MAX_NODES_PER_ZONE][3];
for (int n=0; n<outcell.nnodes; n++)
{
int id = outcell.ids[n];
if (id < mir.origNPoints)
{
coords[n][0] = mir.origXCoords[id];
coords[n][1] = mir.origYCoords[id];
coords[n][2] = mir.origZCoords[id];
}
else
{
coords[n][0] = mir.coordsList[id-mir.origNPoints].x;
coords[n][1] = mir.coordsList[id-mir.origNPoints].y;
coords[n][2] = mir.coordsList[id-mir.origNPoints].z;
}
}
double vol = CalculateVolumeOrAreaHelper(outcell.celltype, coords);
totalvol += vol;
outputvfs[outcell.mat] += vol;
}
}
// And finally, normalize by the total volume
if (outputvfs)
{
for (int matno=0; matno < nMaterials; matno++)
{
outputvfs[matno] /= totalvol;
}
}
}
......@@ -52,6 +52,10 @@
// Programmer: Jeremy Meredith
// Creation: August 18, 2005
//
// Modifications:
// Jeremy Meredith, Fri Feb 13 11:05:06 EST 2009
// Added calculation of output vf's per material, if requested.
//
// ****************************************************************************
class IsovolumeCellReconstructor : public CellReconstructor
......@@ -60,7 +64,7 @@ class IsovolumeCellReconstructor : public CellReconstructor
IsovolumeCellReconstructor(vtkDataSet*, avtMaterial*, ResampledMat&,
int, int, MIRConnectivity&, ZooMIR&);
void ReconstructCell(int, int, int, int*);
void ReconstructCell(int, int, int, int*, double*);
};
#endif
......@@ -111,10 +111,15 @@ RecursiveCellReconstructor::RecursiveCellReconstructor(vtkDataSet *d,
// Mark C. Miller, Thu Feb 9 21:06:10 PST 2006
// Renamed Array class to VisItArray to avoid name collisions with
// third-party libs
//
// Jeremy Meredith, Fri Feb 13 11:05:06 EST 2009
// Added calculation of output vf's per material, if requested.
//
// ****************************************************************************
void
RecursiveCellReconstructor::ReconstructCell(int cellid_, int celltype_,
int nids_, int *ids_)
int nids_, int *ids_,
double *outputvfs)
{
cellid = cellid_;
celltype = celltype_;
......@@ -433,6 +438,17 @@ RecursiveCellReconstructor::ReconstructCell(int cellid_, int celltype_,
}
}
// If we're going to calculate actual volume fractions, first
// zero them out, then accumulate the output cell partial contributions.
double totalvol = 0;
if (outputvfs)
{
for (int matno=0; matno < nMaterials; matno++)
{
outputvfs[matno] = 0.0;
}
}
// Spit the reconstructed cells into the output zone list
for (int out = 0 ; out < outlist.size() ; out++)
{
......@@ -449,6 +465,40 @@ RecursiveCellReconstructor::ReconstructCell(int cellid_, int celltype_,
for (int n=0; n<outcell.nnodes; n++)
mir.indexList.push_back(outcell.ids[n]);
if (outputvfs)
{
double coords[MAX_NODES_PER_ZONE][3];
for (int n=0; n<outcell.nnodes; n++)
{
int id = outcell.ids[n];
if (id < mir.origNPoints)
{
coords[n][0] = mir.origXCoords[id];
coords[n][1] = mir.origYCoords[id];
coords[n][2] = mir.origZCoords[id];
}
else
{
coords[n][0] = mir.coordsList[id-mir.origNPoints].x;
coords[n][1] = mir.coordsList[id-mir.origNPoints].y;
coords[n][2] = mir.coordsList[id-mir.origNPoints].z;
}
}
double vol = CalculateVolumeOrAreaHelper(outcell.celltype, coords);
totalvol += vol;
outputvfs[outcell.mat] += vol;
}
}
// And finally, normalize by the total volume
if (outputvfs)
{
for (int matno=0; matno < nMaterials; matno++)
{
outputvfs[matno] /= totalvol;
}
}
}
......@@ -52,6 +52,10 @@
// Programmer: Jeremy Meredith
// Creation: August 18, 2005
//
// Modifications:
// Jeremy Meredith, Fri Feb 13 11:05:06 EST 2009
// Added calculation of output vf's per material, if requested.
//
// ****************************************************************************
class RecursiveCellReconstructor : public CellReconstructor
......@@ -60,7 +64,7 @@ class RecursiveCellReconstructor : public CellReconstructor
RecursiveCellReconstructor(vtkDataSet*, avtMaterial*, ResampledMat&,
int, int, MIRConnectivity&, ZooMIR&);
void ReconstructCell(int, int, int, int*);
void ReconstructCell(int, int, int, int*, double*);
};
#endif
......@@ -210,9 +210,18 @@ ZooMIR::ReconstructMesh(vtkDataSet *mesh_orig, avtMaterial *mat_orig, int dim)
MIRConnectivity conn;
conn.SetUpConnectivity(mesh);
// extract coordinate arrays
SetUpCoords();
// Pack the material
avtMaterial *mat = mat_orig->CreatePackedMaterial();
// For iteration, create another copy of the material to use
// for the desired VF's. CreatePackedMaterial efficiently
// creates a copy of an already-packed material.
avtMaterial *matCopy = options.numIterations > 0 ?
mat->CreatePackedMaterial() : NULL;
mapMatToUsedMat = mat_orig->GetMapMatToUsedMat();
mapUsedMatToMat = mat_orig->GetMapUsedMatToMat();
......@@ -223,49 +232,81 @@ ZooMIR::ReconstructMesh(vtkDataSet *mesh_orig, avtMaterial *mat_orig, int dim)
int nCells = mesh->GetNumberOfCells();
int nPoints = mesh->GetNumberOfPoints();
// extract coordinate arrays
SetUpCoords();
// resample the material volume fractions to the nodes
ResampledMat rm(nCells, nPoints, mat, &conn);
rm.Resample();
for (int iter = 0; iter < options.numIterations+1; iter++)
{
coordsList.clear();
zonesList.clear();
indexList.clear();
// loop over the cells and do the reconstruction
int timerHandle2 = visitTimer->StartTimer();
ResampledMat rm(nCells, nPoints, mat, &conn);
rm.Resample();
coordsList.reserve(nPoints/4);
zonesList.reserve(int(float(nCells)*1.5));
indexList.reserve(int(float(nCells)*10));
// loop over the cells and do the reconstruction
int timerHandle2 = visitTimer->StartTimer();
CellReconstructor *cr;
if (options.algorithm == 1) // Recursive Clipping
{
cr = new RecursiveCellReconstructor(mesh, mat, rm, nPoints,
nCells, conn, *this);
}
else // algorithm == 2 // Isovolume
{
cr = new IsovolumeCellReconstructor(mesh, mat, rm, nPoints,
nCells, conn, *this);
}
coordsList.reserve(nPoints/4);
zonesList.reserve(int(float(nCells)*1.5));
indexList.reserve(int(float(nCells)*10));
CellReconstructor *cr;
if (options.algorithm == 1) // Recursive Clipping
{
cr = new RecursiveCellReconstructor(mesh, mat, rm, nPoints,
nCells, conn, *this);
}
else // algorithm == 2 // Isovolume
{
cr = new IsovolumeCellReconstructor(mesh, mat, rm, nPoints,
nCells, conn, *this);
}
int *conn_ptr = conn.connectivity;
double actualVFStorage[nMaterials];
double *actualVFs = (options.numIterations>0) ? actualVFStorage : NULL;
double maxdiff = 0;
for (int c = 0 ; c < nCells ; c++, conn_ptr += (*conn_ptr) + 1)
{
int celltype = conn.celltype[c];
int *ids = conn_ptr+1;
int nids = *conn_ptr;
int *conn_ptr = conn.connectivity;
for (int c = 0 ; c < nCells ; c++, conn_ptr += (*conn_ptr) + 1)
{
int celltype = conn.celltype[c];
int *ids = conn_ptr+1;
int nids = *conn_ptr;
cr->ReconstructCell(c, celltype, nids, ids, actualVFs);
cr->ReconstructCell(c, celltype, nids, ids);
}
if (options.numIterations > 0 &&
iter < options.numIterations)
{
std::vector<float> neededVFs, triedVFs;
matCopy->GetVolFracsForZone(c,neededVFs);
mat->GetVolFracsForZone(c,triedVFs);
for (int m=0; m<nMaterials; m++)
{
float wanted = neededVFs[m];
float tried = triedVFs[m];
float got = actualVFs[m];
float diff = wanted - got;
if (diff != 0)
{
float newval = tried + options.iterationDamping * diff;
mat->SetVolFracForZoneAndMat(c,m,newval);
if (diff > maxdiff)
maxdiff = diff;
}
}
}
}