Updates will be applied April 15th at 12pm EDT (UTC-0400). GitLab could be a little slow between 12 - 12:45pm EDT.

Commit 5fd700cc authored by hrchilds's avatar hrchilds

Merge John Anderson's material interface reconstruction algorithm into

the trunk.  Algorithm is called Discrete.


git-svn-id: http://visit.ilight.com/svn/visit/trunk/src@9052 18c085ea-50e0-402c-830e-de6fd14e8384
parent 2b79ded8
......@@ -99,6 +99,7 @@
#include <PickAttributes.h>
#include <PickVarInfo.h>
#ifndef DBIO_ONLY
#include <DiscreteMIR.h>
#include <TetMIR.h>
#include <ZooMIR.h>
#include <YoungsMIR.h>
......@@ -3476,11 +3477,15 @@ avtGenericDatabase::AddOriginalNodesArray(vtkDataSet *ds, const int domain)
// they had only one material. Added code to account for that in
// case of creating boundary surfaces.
//
// John C. Anderson, Thu Jan 15 10:20:20 2009
// Added annealing time for Discrete MIR.
//
// Jeremy Meredith, Fri Feb 13 11:22:39 EST 2009
// Added MIR iteration capability.
//
// Mark C. Miller, Wed Mar 4 18:00:22 PST 2009
// Adjusted for dbio-only build
//
// ****************************************************************************
avtDataTree_p
......@@ -3499,6 +3504,7 @@ avtGenericDatabase::MaterialSelect(vtkDataSet *ds, avtMaterial *mat,
int mirNumIterations,
float mirIterationDamping,
float isovolumeMIRVF,
int annealingTime,
bool didGhosts,
bool &subdivisionOccurred,
bool &notAllCellsSubdivided,
......@@ -3564,6 +3570,7 @@ avtGenericDatabase::MaterialSelect(vtkDataSet *ds, avtMaterial *mat,
mirNumIterations,
mirIterationDamping,
isovolumeMIRVF,
annealingTime,
didGhosts,
subdivisionOccurred,
notAllCellsSubdivided,
......@@ -4363,6 +4370,12 @@ avtGenericDatabase::SpeciesSelect(avtDatasetCollection &dsc,
// Do not delete the new material if we do "simplify heavily mixed",
// since it might be used later.
//
// John C. Anderson, Fri Oct 17 17:37:16 2008
// Added Discrete MIR option.
//
// John C. Anderson, Thu Jan 15 10:20:20 2009
// Added annealing time for Discrete MIR.
//
// Hank Childs, Thu Feb 21 16:41:25 PST 2008
// Throw an exception if we get a bad MIR type value.
//
......@@ -4394,6 +4407,7 @@ avtGenericDatabase::GetMIR(int domain, const char *varname, int timestep,
int mirNumIterations,
float mirIterationDamping,
float isovolumeMIRVF,
int annealingTime,
bool didGhosts,
bool &subdivisionOccurred,
bool &notAllCellsSubdivided, bool reUseMIR,
......@@ -4417,7 +4431,7 @@ avtGenericDatabase::GetMIR(int domain, const char *varname, int timestep,
}
char cacheLbl[1000];
sprintf(cacheLbl, "MIR_%s_%s_%s_%s_%s_%d_%f_%s_%d_%f",
sprintf(cacheLbl, "MIR_%s_%s_%s_%s_%s_%d_%f_%s_%d_%d_%f",
needValidConnectivity ? "FullSubdiv" : "MinimalSubdiv",
needSmoothMaterialInterfaces ? "Smooth" : "NotSmooth",
needCleanZonesOnly ? "CleanOnly" : "SplitMixed",
......@@ -4428,11 +4442,12 @@ avtGenericDatabase::GetMIR(int domain, const char *varname, int timestep,
mirAlgorithm==0 ? "TetMIR" :
(mirAlgorithm==1 ? "ZooMIR" :
(mirAlgorithm==2 ? "IsovolumeMIR" :
"YoungsMIR")),
(mirAlgorithm==3 ? "YoungsMIR" :
"DiscreteMIR"))),
annealingTime,
mirNumIterations,
mirIterationDamping);
//
// See if we already have the data lying around.
//
......@@ -4476,6 +4491,10 @@ avtGenericDatabase::GetMIR(int domain, const char *varname, int timestep,
mir = new YoungsMIR;
break;
case 4:
mir = new DiscreteMIR;
break;
default:
EXCEPTION0(ImproperUseException);
}
......@@ -4487,6 +4506,7 @@ avtGenericDatabase::GetMIR(int domain, const char *varname, int timestep,
mir->SetSmoothing(needSmoothMaterialInterfaces);
mir->SetCleanZonesOnly(needCleanZonesOnly);
mir->SetIsovolumeVF(isovolumeMIRVF);
mir->SetAnnealingTime(annealingTime);
if (topoDim == 3)
{
mir->Reconstruct3DMesh(ds, mat_to_use);
......@@ -7376,6 +7396,9 @@ avtGenericDatabase::ApplyGhostForDomainNesting(avtDatasetCollection &ds,
// Hank Childs, Wed Jul 25 14:16:36 PDT 2007
// Renamed method: NeedBoundarySurfaces -> GetBoundarySurfaceRepresentation.
//
// John C. Anderson, Thu Jan 15 10:20:20 2009
// Added annealing time for Discrete MIR.
//
// Jeremy Meredith, Fri Feb 13 11:22:39 EST 2009
// Added MIR iteration capability.
//
......@@ -7456,6 +7479,7 @@ avtGenericDatabase::MaterialSelect(avtDatasetCollection &ds,
spec->MIRNumIterations(),
spec->MIRIterationDamping(),
spec->IsovolumeMIRVF(),
spec->AnnealingTime(),
didGhosts, so, nacs, reUseMIR);
notAllCellsSubdivided = notAllCellsSubdivided || nacs ||
......
......@@ -322,6 +322,9 @@ class vtkUnstructuredGrid;
// Hank Childs, Sun Oct 28 21:09:44 PST 2007
// Added Boolean argument to GetDomainBoundaryInformation.
//
// John C. Anderson, Thu Jan 15 10:20:20 2009
// Added annealing time to GetMIR and MaterialSelect.
//
// Hank Childs, Sun Feb 10 19:43:59 MST 2008
// Add support for streaming ghost generation.
//
......@@ -454,12 +457,12 @@ class DATABASE_API avtGenericDatabase : public avtDatasetDatabase
stringVector &,
bool, bool, bool, bool, bool,
bool, bool, int, int,
int, float, float, bool,
int, float, float, int, bool,
bool&, bool&, bool);
void_ref_ptr GetMIR(int, const char *, int, vtkDataSet*,
avtMaterial *, int, bool, bool, bool,
bool, int, int, int, float,
float, bool, bool&, bool&,bool,
float, int, bool, bool&, bool&,bool,
avtMaterial *&);
avtMaterial *GetMaterial(int, const char *, int, const avtDataRequest_p = 0);
avtSpecies *GetSpecies(int, const char *, int);
......
......@@ -199,6 +199,22 @@ MIR::SetIsovolumeVF(float vf)
options.isovolumeVF = vf;
}
// ****************************************************************************
// Method: MIR::SetAnnealingTime
//
// Purpose:
// Set the maximum simulated annealing time for Discrete MIR.
//
// Programmer: John C. Anderson
// Creation: January 15, 2009
//
// ****************************************************************************
void
MIR::SetAnnealingTime(int t)
{
options.annealingTime = t;
}
// ****************************************************************************
// Default Constructor: MIR::MIR
//
......
......@@ -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.
//
// John C. Anderson, Thu Jan 15 10:20:20 2009
// Added annealing time for Discrete MIR.
//
// Jeremy Meredith, Fri Feb 13 11:22:39 EST 2009
// Added MIR iteration capability.
//
......@@ -99,6 +102,7 @@ class MIR_API MIR
void SetLeaveCleanZonesWhole(bool);
void SetCleanZonesOnly(bool);
void SetIsovolumeVF(float);
void SetAnnealingTime(int);
// do the processing
virtual bool Reconstruct3DMesh(vtkDataSet *, avtMaterial *) = 0;
......
......@@ -48,6 +48,9 @@
// Jeremy Meredith, Thu Aug 18 16:36:55 PDT 2005
// Added algorithm and isovolumeVF.
//
// John C. Anderson, Thu Jan 15 10:20:20 2009
// Added annealing time for Discrete MIR.
//
// Jeremy Meredith, Fri Feb 13 11:22:39 EST 2009
// Added MIR iteration capability.
//
......@@ -62,4 +65,5 @@ MIROptions::MIROptions()
leaveCleanZonesWhole = true;
cleanZonesOnly = false;
isovolumeVF = 0.5;
annealingTime = 10;
}
......@@ -75,6 +75,9 @@
// Jeremy Meredith, Thu Aug 18 16:36:42 PDT 2005
// Added algorithm and isovolumeVF.
//
// John C. Anderson, Thu Jan 15 10:20:20 2009
// Added annealing time for Discrete MIR.
//
// Jeremy Meredith, Fri Feb 13 11:22:39 EST 2009
// Added MIR iteration and damping capability.
//
......@@ -97,6 +100,7 @@ class MIR_API MIROptions
bool leaveCleanZonesWhole;
bool cleanZonesOnly;
float isovolumeVF;
int annealingTime;
public:
MIROptions();
......
This diff is collapsed.
/*****************************************************************************
*
* Copyright (c) 2000 - 2008, Lawrence Livermore National Security, LLC
* Produced at the Lawrence Livermore National Laboratory
* LLNL-CODE-400142
* All rights reserved.
*
* This file is part of VisIt. For details, see https://visit.llnl.gov/. The
* full copyright notice is contained in the file COPYRIGHT located at the root
* of the VisIt distribution or at http://www.llnl.gov/visit/copyright.html.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
*
* - Redistributions of source code must retain the above copyright notice,
* this list of conditions and the disclaimer below.
* - Redistributions in binary form must reproduce the above copyright notice,
* this list of conditions and the disclaimer (as noted below) in the
* documentation and/or other materials provided with the distribution.
* - Neither the name of the LLNS/LLNL nor the names of its contributors may
* be used to endorse or promote products derived from this software without
* specific prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL LAWRENCE LIVERMORE NATIONAL SECURITY,
* LLC, THE U.S. DEPARTMENT OF ENERGY OR CONTRIBUTORS BE LIABLE FOR ANY
* DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
* SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
* CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
* DAMAGE.
*
*****************************************************************************/
#ifndef DISCRETEMIR_H
#define DISCRETEMIR_H
#include <MIR.h>
#include <MIRConnectivity.h>
#include <VisItArray.h>
// ****************************************************************************
// Class: DiscreteMIR
//
// Purpose:
// Produces a discretized material reconstruction based upon:
//
// John C. Anderson, Christoph Garth, Mark A. Duchaineau, and
// Kenneth I. Joy, "Discrete Multi-Material Interface
// Reconstruction for Volume Fraction Data," Computer Graphics
// Forum (Proc. of EuroVis), vol. 27, pp. 1015--1022, May 2008.
//
// Notes:
//
// I haven't looked very hard, but I could probably eliminate the
// Cell and Node classes used here in favor of code already
// included withing VTK or VisIt. For now, however, they work...
//
// Programmer: John C. Anderson
// Creation: October 17, 2008
//
// Modifications:
//
// ****************************************************************************
class MIR_API DiscreteMIR : public MIR
{
public:
DiscreteMIR();
virtual ~DiscreteMIR();
// do the processing
bool ReconstructMesh(vtkDataSet *, avtMaterial *, int);
virtual bool Reconstruct3DMesh(vtkDataSet *, avtMaterial *);
virtual bool Reconstruct2DMesh(vtkDataSet *, avtMaterial *);
// material select everything -- all variables, the mesh, and the material
// if requested.
virtual vtkDataSet *GetDataset(std::vector<int>, vtkDataSet *,
std::vector<avtMixedVariable *>, bool,
avtMaterial * = NULL);
// get some result flags
virtual bool SubdivisionOccurred() { return true; }
virtual bool NotAllCellsSubdivided() { return false; }
protected:
struct ReconstructedCoord
{
double x,y,z;
double weight[MAX_NODES_PER_ZONE];
int origzone;
};
struct ReconstructedZone
{
int origzone;
int startindex;
int mix_index;
int mat;
unsigned char celltype;
unsigned char nnodes;
};
int origNPoints;
std::vector<float> origXCoords;
std::vector<float> origYCoords;
std::vector<float> origZCoords;
VisItArray<ReconstructedCoord> coordsList;
VisItArray<ReconstructedZone> zonesList;
VisItArray<vtkIdType> indexList;
int dimension;
int nMaterials;
int nOrigMaterials;
bool noMixedZones;
int singleMat;
std::vector<int> mapMatToUsedMat;
std::vector<int> mapUsedMatToMat;
vtkPoints *outPts;
vtkDataSet *mesh;
protected:
bool ReconstructCleanMesh(vtkDataSet *, avtMaterial *);
void SetUpCoords();
protected: // John's stuff...
class Node;
class Cell
{
public:
Cell():
m_i(-1),
m_j(-1),
m_k(0)
{
}
Cell(int _m_i, int _m_j, int _m_k = 0):
m_i(_m_i),
m_j(_m_j),
m_k(_m_k)
{
}
#if 0
inline seagull::Point3 centroid4() const
{
seagull::Point3 centroid(0.0f, 0.0f, 0.0f);
for(int n = 0; n < 4; ++n)
{
Gallimaufry::Node node = incidentNode4(n);
centroid += seagull::Point3(node.m_i, node.m_j, node.m_k);
}
centroid /= 4;
return centroid;
}
inline seagull::Point3 centroid8() const
{
seagull::Point3 centroid(0.0f, 0.0f, 0.0f);
for(int n = 0; n < 8; ++n)
{
Gallimaufry::Node node = incidentNode8(n);
centroid += seagull::Point3(node.m_i, node.m_j, node.m_k);
}
centroid /= 8;
return centroid;
}
#endif
inline Cell incidentCell4(int c) const
{
return Cell(m_i + m_incidentCell4[c][0],
m_j + m_incidentCell4[c][1]);
}
inline Cell incidentCell6(int c) const
{
return Cell(m_i + m_incidentCell6[c][0],
m_j + m_incidentCell6[c][1],
m_k + m_incidentCell6[c][2]);
}
inline Cell incidentCell8(int c) const
{
return Cell(m_i + m_incidentCell8[c][0],
m_j + m_incidentCell8[c][1]);
}
inline Cell incidentCell26(int c) const
{
return Cell(m_i + m_incidentCell26[c][0],
m_j + m_incidentCell26[c][1],
m_k + m_incidentCell26[c][2]);
}
inline Node incidentNode4(int n) const
{
return Node(m_i + m_incidentNode4[n][0],
m_j + m_incidentNode4[n][1]);
}
inline Node incidentNode8(int n) const
{
return Node(m_i + m_incidentNode8[n][0],
m_j + m_incidentNode8[n][1],
m_k + m_incidentNode8[n][2]);
}
bool operator==(const Cell &that) const
{
return (m_i == that.m_i &&
m_j == that.m_j &&
m_k == that.m_k);
}
bool operator!=(const Cell &that) const
{
return !operator==(that);
}
bool operator<(const Cell &that) const
{
if(m_i < that.m_i)
return true;
else if(m_i == that.m_i &&
m_j < that.m_j)
return true;
else if(m_i == that.m_i &&
m_j == that.m_j &&
m_k < that.m_k)
return true;
else
return false;
}
size_t operator()(const Cell &that) const
{
return that.m_i * that.m_j * that.m_k;
}
int m_i, m_j, m_k;
protected:
// Clockwise enumeration of 4 face-incident cells.
static int m_incidentCell4[4][2];
/// Face-incident cells in 3D.
static int m_incidentCell6[6][3];
/// Clockwise enumeration of 8 incident cells.
static int m_incidentCell8[8][2];
/// Incident cells in 3D.
static int m_incidentCell26[26][3];
/// Incident nodes in 2D.
static int m_incidentNode4[4][2];
/// Incident nodes in 3D.
static int m_incidentNode8[8][3];
};
class Node
{
public:
Node():
m_i(-1),
m_j(-1),
m_k(0)
{
}
Node(int _m_i, int _m_j, int _m_k = 0):
m_i(_m_i),
m_j(_m_j),
m_k(_m_k)
{
}
inline Cell incidentCell4(int c) const
{
return Cell(m_i + m_incidentCell4[c][0],
m_j + m_incidentCell4[c][1]);
}
inline Cell incidentCell8(int c) const
{
return Cell(m_i + m_incidentCell8[c][0],
m_j + m_incidentCell8[c][1],
m_k + m_incidentCell8[c][2]);
}
inline Node neighborNode(int n) const
{
return Node(m_i + m_neighborNode[n][0],
m_j + m_neighborNode[n][1]);
}
bool operator==(const Node &that) const
{
return (m_i == that.m_i &&
m_j == that.m_j &&
m_k == that.m_k);
}
bool operator!=(const Node &that) const
{
return !operator==(that);
}
bool operator<(const Node &that) const
{
if(m_i < that.m_i)
return true;
else if(m_i == that.m_i &&
m_j < that.m_j)
return true;
else if(m_i == that.m_i &&
m_j == that.m_j &&
m_k < that.m_k)
return true;
else
return false;
}
size_t operator()(const Node &that) const
{
return that.m_i * that.m_j * that.m_k;
}
int m_i, m_j, m_k;
protected:
static int m_incidentCell4[4][2];
static int m_incidentCell8[8][3];
static int m_neighborNode[4][2];
};
int DX, DY, DZ;
int NX, NY, NZ;
int dimensions[3];
float *xspacing, *yspacing, *zspacing;
std::vector< Cell > m_mixedCells;
unsigned char **m_labels;
float *m_neighborhood;
size_t m_size[3];
double m_temperature;
void optimize();
int id(const Cell &cell) const
{
return
cell.m_k * (dimensions[0] * dimensions[1]) +
cell.m_j * (dimensions[0]) +
cell.m_i;
}
int id(const Node &node) const
{
return
node.m_k * ((dimensions[0] + 1) * (dimensions[1] + 1)) +
node.m_j * ((dimensions[0] + 1)) +
node.m_i;
}
bool isValid(const Cell &cell) const;
bool isValid(const Node &node) const;
size_t getSize(int d) const
{
return m_size[d];
}
unsigned char get(size_t i, size_t j, size_t k) const;
bool inside(size_t i, size_t j, size_t k) const
{
return
i