Commit a4cc7836 authored by Philippe Pébaÿ's avatar Philippe Pébaÿ
Browse files

Fixed issues with the computation of material interfaces in HTG

This commit addresses to main issues:
1. The AxisReflection was /adding/ modified arrays, as opposed to
   modifying in situ the interface arrays (vectors and intercepts).
   This had initially made on purpose with the intent of temporary
   retaining the original arrays for debugging purposes. However,
   this had been forgotten and left in place. As a result, interface
   computations after a reflection where incorrect. This bug was
   not revealing itself in the test harness so far as this combination
   (reflection + interface computations) was not exercised.
   As a result, this commit addresses the issue but also verifies
   this combination. Furthermore, the HTG source was modified in
   order to allow for the specification of intercepts arrays in
   the form that exercise this configuration.
2. The Geometry extraction was not building material interfaces, as
   a result of some code having been leftover during the HTG-v2
   merger in VTK 8. This commit fixes this issue, and also adds
   a unit test to verify material extraction together with its
   corresponding non-regression image.
parent 60a2e996
......@@ -720,8 +720,8 @@ protected:
bool InitPureMaterialMask;
bool HasInterface;
char *InterfaceNormalsName;
char *InterfaceInterceptsName;
char* InterfaceNormalsName;
char* InterfaceInterceptsName;
vtkDataArray* XCoordinates;
vtkDataArray* YCoordinates;
......
......@@ -20,6 +20,7 @@ vtk_add_test_cxx(${vtk-module}CxxTests tests
TestHyperTreeGridBinary2DContourMaterial.cxx
TestHyperTreeGridBinary2DDepthLimiter.cxx
TestHyperTreeGridBinary2DDepthLimiterMaterial.cxx
TestHyperTreeGridBinary2DInterfaceMaterial.cxx
TestHyperTreeGridBinary2DMaterial.cxx
TestHyperTreeGridBinary2DMaterialIJK.cxx
TestHyperTreeGridBinary2DThreshold.cxx
......
/*==================================================================
Program: Visualization Toolkit
Module: TestHyperTreeGridBinary2DInterfaceMaterial.cxx
Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
All rights reserved.
See Copyright.txt or http://www.kitware.com/Copyright.htm 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.
===================================================================*/
// .SECTION Thanks
// This test was written by Philippe Pebay, NexGen Analytics 2017
// This work was supported by Commissariat a l'Energie Atomique (CEA/DIF)
#include "vtkHyperTreeGrid.h"
#include "vtkHyperTreeGridGeometry.h"
#include "vtkHyperTreeGridSource.h"
#include "vtkCamera.h"
#include "vtkCellData.h"
#include "vtkNew.h"
#include "vtkPointData.h"
#include "vtkPolyData.h"
#include "vtkPolyDataMapper.h"
#include "vtkProperty.h"
#include "vtkRegressionTestImage.h"
#include "vtkRenderer.h"
#include "vtkRenderWindow.h"
#include "vtkRenderWindowInteractor.h"
int TestHyperTreeGridBinary2DInterfaceMaterial( int argc, char* argv[] )
{
// Hyper tree grid
vtkNew<vtkHyperTreeGridSource> htGrid;
htGrid->SetMaximumLevel( 6 );
htGrid->SetDimension( 2 );
htGrid->SetOrientation( 2 ); // in xy plane
htGrid->SetGridSize( 2, 3, 1 );
htGrid->SetGridScale( 1.5, 1., 10. ); // this is to test that orientation fixes scale
htGrid->SetBranchFactor( 2 );
htGrid->SetDescriptor( "RRRRR.|.... .R.. RRRR R... R...|.R.. ...R ..RR .R.. R... .... ....|.... ...R ..R. .... .R.. R...|.... .... .R.. ....|...." );
htGrid->GenerateInterfaceFieldsOn();
htGrid->Update();
vtkHyperTreeGrid* H = vtkHyperTreeGrid::SafeDownCast( htGrid->GetOutput() );
H->SetHasInterface( 1 );
char normalsName[] = "Normals";
H->SetInterfaceNormalsName( normalsName );
char interceptsName[] = "Intercepts";
H->SetInterfaceInterceptsName( interceptsName );
// Modify intercepts array
vtkDataArray* interArray = vtkDataSet::SafeDownCast( htGrid->GetOutput() )->GetPointData()->GetArray( "Intercepts" );
for ( vtkIdType i = 0; i < interArray->GetNumberOfTuples(); ++ i )
{
double* inter = interArray->GetTuple3( i );
interArray->SetTuple3( i, -.25, -.5, -1. );
}
// Geometries
vtkNew<vtkHyperTreeGridGeometry> geometry1;
geometry1->SetInputConnection( htGrid->GetOutputPort() );
geometry1->Update();
vtkPolyData* pd = geometry1->GetPolyDataOutput();
vtkNew<vtkHyperTreeGridGeometry> geometry2;
geometry2->SetInputConnection( htGrid->GetOutputPort() );
// Mappers
vtkMapper::SetResolveCoincidentTopologyToPolygonOffset();
vtkNew<vtkPolyDataMapper> mapper1;
mapper1->SetInputConnection( geometry1->GetOutputPort() );
mapper1->ScalarVisibilityOff();
vtkNew<vtkPolyDataMapper> mapper2;
mapper2->SetInputConnection( geometry2->GetOutputPort() );
mapper2->SetScalarRange( pd->GetCellData()->GetScalars()->GetRange() );
// Actors
vtkNew<vtkActor> actor1;
actor1->SetMapper( mapper1 );
actor1->GetProperty()->SetRepresentationToWireframe();
actor1->GetProperty()->SetColor( .7, .7, .7 );
vtkNew<vtkActor> actor2;
actor2->SetMapper( mapper2 );
// Camera
double bd[6];
pd->GetBounds( bd );
vtkNew<vtkCamera> camera;
camera->SetClippingRange( 1., 100. );
camera->SetFocalPoint( pd->GetCenter() );
camera->SetPosition( .5 * bd[1], .5 * bd[3], 6. );
// Renderer
vtkNew<vtkRenderer> renderer;
renderer->SetActiveCamera( camera );
renderer->SetBackground( 1., 1., 1. );
renderer->AddActor( actor1 );
renderer->AddActor( actor2 );
// Render window
vtkNew<vtkRenderWindow> renWin;
renWin->AddRenderer( renderer );
renWin->SetSize( 400, 400 );
renWin->SetMultiSamples( 0 );
// Interactor
vtkNew<vtkRenderWindowInteractor> iren;
iren->SetRenderWindow( renWin );
// Render and test
renWin->Render();
int retVal = vtkRegressionTestImageThreshold( renWin, 70 );
if ( retVal == vtkRegressionTester::DO_INTERACTOR )
{
iren->Start();
}
return !retVal;
}
......@@ -14,6 +14,7 @@
===================================================================*/
// .SECTION Thanks
// This test was written by Philippe Pebay, 2016
// This test was modified by Philippe Pebay, NexGen Analytics 2017
// This work was supported by Commissariat a l'Energie Atomique (CEA/DIF)
#include "vtkHyperTreeGrid.h"
......@@ -38,21 +39,20 @@ int TestHyperTreeGridBinary2DVector( int argc, char* argv[] )
{
// Hyper tree grid
vtkNew<vtkHyperTreeGridSource> htGrid;
int maxLevel = 6;
htGrid->SetMaximumLevel( maxLevel );
htGrid->SetMaximumLevel( 6 );
htGrid->SetOrientation( 2 ); // in xy plane
htGrid->SetBranchFactor( 2 );
htGrid->SetGridSize( 2, 3, 1 );
htGrid->SetGridScale( 1.5, 1., 10. ); // this is to test that orientation fixes scale
htGrid->SetDimension( 2 );
htGrid->SetOrientation( 2 ); // in xy plane
htGrid->SetBranchFactor( 2 );
htGrid->SetDescriptor( "RRRRR.|.... .R.. RRRR R... R...|.R.. ...R ..RR .R.. R... .... ....|.... ...R ..R. .... .R.. R...|.... .... .R.. ....|...." );
htGrid->GenerateVectorFieldOn();
htGrid->GenerateInterfaceFieldsOn();
htGrid->Update();
vtkHyperTreeGrid* H = vtkHyperTreeGrid::SafeDownCast( htGrid->GetOutput() );
H->SetHasInterface( 1 );
char normalsName[] = "Vector";
char normalsName[] = "Normals";
H->SetInterfaceNormalsName( normalsName );
char interceptsName[] = "Depth";
char interceptsName[] = "Intercepts";
H->SetInterfaceInterceptsName( interceptsName );
// Cell centers
......
......@@ -14,6 +14,7 @@
===================================================================*/
// .SECTION Thanks
// This test was written by Philippe Pebay, 2016
// This test was modified by Philippe Pebay, NexGen Analytics 2017
// This work was supported by Commissariat a l'Energie Atomique (CEA/DIF)
#include "vtkHyperTreeGrid.h"
......@@ -39,20 +40,21 @@ int TestHyperTreeGridBinary2DVectorAxisReflectionXCenter( int argc, char* argv[]
{
// Hyper tree grid
vtkNew<vtkHyperTreeGridSource> htGrid;
int maxLevel = 6;
htGrid->SetMaximumLevel( maxLevel );
htGrid->SetGridSize( 2, 3, 1 );
htGrid->SetGridScale( 1.5, 1., 10. ); // this is to test that orientation fixes scale
htGrid->SetMaximumLevel( 6 );
htGrid->SetDimension( 2 );
htGrid->SetOrientation( 2 ); // in xy plane
htGrid->SetGridSize( 2, 3, 1 );
htGrid->SetGridScale( 1.5, 1., 10. ); // this is to test that orientation fixes scale
htGrid->SetBranchFactor( 2 );
htGrid->SetDescriptor( "RRRRR.|.... .R.. RRRR R... R...|.R.. ...R ..RR .R.. R... .... ....|.... ...R ..R. .... .R.. R...|.... .... .R.. ....|...." );
htGrid->GenerateVectorFieldOn();
htGrid->GenerateInterfaceFieldsOn();
htGrid->Update();
vtkHyperTreeGrid* H = vtkHyperTreeGrid::SafeDownCast( htGrid->GetOutput() );
H->SetHasInterface( 1 );
H->SetInterfaceNormalsName( "Vector" );
H->SetInterfaceInterceptsName( "Depth" );
char normalsName[] = "Normals";
H->SetInterfaceNormalsName( normalsName );
char interceptsName[] = "Intercepts";
H->SetInterfaceInterceptsName( interceptsName );
// Axis reflection
vtkNew<vtkHyperTreeGridAxisReflection> reflection;
......
......@@ -14,6 +14,7 @@
===================================================================*/
// .SECTION Thanks
// This test was written by Philippe Pebay, 2016
// This test was modified by Philippe Pebay, NexGen Analytics 2017
// This work was supported by Commissariat a l'Energie Atomique (CEA/DIF)
#include "vtkHyperTreeGrid.h"
......@@ -39,21 +40,20 @@ int TestHyperTreeGridBinary2DVectorAxisReflectionYCenter( int argc, char* argv[]
{
// Hyper tree grid
vtkNew<vtkHyperTreeGridSource> htGrid;
int maxLevel = 6;
htGrid->SetMaximumLevel( maxLevel );
htGrid->SetGridSize( 2, 3, 1 );
htGrid->SetGridScale( 1.5, 1., 10. ); // this is to test that orientation fixes scale
htGrid->SetMaximumLevel( 6 );
htGrid->SetDimension( 2 );
htGrid->SetOrientation( 2 ); // in xy plane
htGrid->SetGridSize( 2, 3, 1 );
htGrid->SetGridScale( 1.5, 1., 10. ); // this is to test that orientation fixes scale
htGrid->SetBranchFactor( 2 );
htGrid->SetDescriptor( "RRRRR.|.... .R.. RRRR R... R...|.R.. ...R ..RR .R.. R... .... ....|.... ...R ..R. .... .R.. R...|.... .... .R.. ....|...." );
htGrid->GenerateVectorFieldOn();
htGrid->GenerateInterfaceFieldsOn();
htGrid->Update();
vtkHyperTreeGrid* H = vtkHyperTreeGrid::SafeDownCast( htGrid->GetOutput() );
H->SetHasInterface( 1 );
char normalsName[] = "Vector";
char normalsName[] = "Normals";
H->SetInterfaceNormalsName( normalsName );
char interceptsName[] = "Depth";
char interceptsName[] = "Intercepts";
H->SetInterfaceInterceptsName( interceptsName );
// Axis reflection
......
......@@ -144,19 +144,6 @@ int vtkHyperTreeGridAxisReflection::ProcessTrees( vtkHyperTreeGrid* input,
vtkDoubleArray* outCoords = vtkDoubleArray::New();
outCoords->SetNumberOfTuples( size );
// Create arrays for reflected interface if present
vtkDoubleArray* outNormals = 0;
vtkDoubleArray* outIntercepts = 0;
if ( hasInterface )
{
vtkIdType nTuples = inNormals->GetNumberOfTuples();
outNormals = vtkDoubleArray::New();
outNormals->SetNumberOfComponents( 3 );
outNormals->SetNumberOfTuples( nTuples );
outIntercepts = vtkDoubleArray::New();
outIntercepts->SetNumberOfTuples( nTuples );
}
// Reflect point coordinate
double coord;
for ( unsigned int i = 0; i < size; ++ i )
......@@ -186,29 +173,19 @@ int vtkHyperTreeGridAxisReflection::ProcessTrees( vtkHyperTreeGrid* input,
for ( vtkIdType i = 0; i < nTuples; ++ i )
{
// Compute and stored reflected normal
double norm[3];
memcpy( norm, inNormals->GetTuple3( i ) , 3 * sizeof( double ) );
double* norm = inNormals->GetTuple3( i );
norm[direction] = - norm[direction];
outNormals->SetTuple3( i, norm[0], norm[1], norm[2] );
inNormals->SetTuple3( i, norm[0], norm[1], norm[2] );
// Compute and store reflected intercept
double inter = inIntercepts->GetTuple1( i );
inter -= 2. * offset * norm[direction];
outIntercepts->SetTuple1( i, inter );
double* inter = inIntercepts->GetTuple3( i );
inter[0] -= 2. * offset * norm[direction];
inIntercepts->SetTuple3( i, inter[0], inter[1], inter[2] );
} // i
// Assign new interface arrays if available
this->OutData->SetVectors( outNormals );
this->OutData->AddArray( outIntercepts );
} // if ( hasInterface )
// Clean up
outCoords->Delete();
if ( hasInterface )
{
outNormals->Delete();
outIntercepts->Delete();
}
return 1;
}
......@@ -24,7 +24,8 @@
* vtkHyperTreeGrid vtkHyperTreeGridAlgorithm vtkReflectionFilter
*
* @par Thanks:
* This class was written by Philippe Pebay on a idea of Guénolé Harel and Jacques-Bernard Lekien, 2016
* This class was written by Philippe Pebay based on a idea of Guenole
* Harel and Jacques-Bernard Lekien, 2016
* This work was supported by Commissariat a l'Energie Atomique (CEA/DIF)
*/
......@@ -115,5 +116,4 @@ private:
void operator=(const vtkHyperTreeGridAxisReflection&) = delete;
};
#endif /* vtkHyperTreeGridAxisReflection */
......@@ -20,9 +20,10 @@
* vtkHyperTreeGrid vtkHyperTreeGridAlgorithm
*
* @par Thanks:
* This class was written by Philippe Pebay, Joachim Pouderoux, and Charles Law, Kitware 2013
* This class was modified by Guénolé Harel and Jacques-Bernard Lekien, 2014
* This class was rewritten by Philippe Pebay, 2016
* This class was written by Philippe Pebay, Joachim Pouderoux,
* and Charles Law, Kitware 2013
* This class was modified by Guenole Harel and Jacques-Bernard Lekien, 2014
* This class was rewritten by Philippe Pebay, NexGen Analytics 2017
* This work was supported by Commissariat a l'Energie Atomique (CEA/DIF)
*/
......@@ -34,8 +35,11 @@
class vtkBitArray;
class vtkCellArray;
class vtkDoubleArray;
class vtkHyperTreeGrid;
class vtkHyperTreeGridCursor;
class vtkIdList;
class vtkIdTypeArray;
class vtkPoints;
class VTKFILTERSHYPERTREE_EXPORT vtkHyperTreeGridGeometry : public vtkHyperTreeGridAlgorithm
......@@ -82,7 +86,7 @@ protected:
/**
* Helper method to generate a face based on its normal and offset from cursor origin
*/
void AddFace( vtkIdType, double*, double*, int, unsigned int );
void AddFace( vtkIdType, double*, double*, int, unsigned int, bool create = true );
/**
* Dimension of input grid
......@@ -104,6 +108,47 @@ protected:
*/
vtkCellArray* Cells;
//@{
/**
* Keep track of input interface parameters
*/
bool HasInterface;
vtkDoubleArray* Normals;
vtkDoubleArray* Intercepts;
//@}
//@{
/**
* Storage for interface points
*/
vtkPoints* FacePoints;
vtkIdList* FaceIDs;
//@}
//@{
/**
* Storage for interface edges
*/
vtkIdType EdgesA[12];
vtkIdType EdgesB[12];
//@}
//@{
/**
* Storage for interface faces
*/
vtkIdTypeArray* FacesA;
vtkIdTypeArray* FacesB;
//@}
//@{
/**
* Storage for interface scalars
*/
vtkDoubleArray* FaceScalarsA;
vtkDoubleArray* FaceScalarsB;
//@}
private:
vtkHyperTreeGridGeometry(const vtkHyperTreeGridGeometry&) = delete;
void operator=(const vtkHyperTreeGridGeometry&) = delete;
......
......@@ -83,8 +83,8 @@ vtkHyperTreeGridSource::vtkHyperTreeGridSource()
// By default do not use the material mask
this->UseMaterialMask = false;
// By default do not generate a vector field
this->GenerateVectorField = false;
// By default do not generate interface vector fields
this->GenerateInterfaceFields = false;
// Grid description & material mask as strings
this->Descriptor = new char[2];
......@@ -204,7 +204,7 @@ void vtkHyperTreeGridSource::PrintSelf( ostream& os, vtkIndent indent )
os << indent << "UseDescriptor: " << this->UseDescriptor << endl;
os << indent << "UseMaterialMask: " << this->UseMaterialMask << endl;
os << indent << "GenerateVectorField:" << this->GenerateVectorField << endl;
os << indent << "GenerateInterfaceFields:" << this->GenerateInterfaceFields << endl;
os << indent << "LevelZeroMaterialIndex: " << this->LevelZeroMaterialIndex << endl;
os << indent << "Descriptor: " << this->Descriptor << endl;
......@@ -527,13 +527,18 @@ int vtkHyperTreeGridSource::RequestData( vtkInformation*,
depthArray->SetNumberOfComponents( 1 );
outData->SetScalars( depthArray );
if ( this->GenerateVectorField )
if ( this->GenerateInterfaceFields )
{
// Prepare array of triples of doubles for vector field
vtkNew<vtkDoubleArray> vectorArray;
vectorArray->SetName( "Vector" );
vectorArray->SetNumberOfComponents( 3 );
outData->SetVectors( vectorArray );
// Prepare arrays of triples for interface surrogates
vtkNew<vtkDoubleArray> normalsArray;
normalsArray->SetName( "Normals" );
normalsArray->SetNumberOfComponents( 3 );
outData->SetVectors( normalsArray );
vtkNew<vtkDoubleArray> interceptsArray;
interceptsArray->SetName( "Intercepts" );
interceptsArray->SetNumberOfComponents( 3 );
outData->AddArray( interceptsArray );
}
if ( ! this->UseDescriptor )
......@@ -864,11 +869,12 @@ void vtkHyperTreeGridSource::SubdivideFromStringDescriptor( vtkHyperTreeGrid* ou
// Set depth array value
outData->GetArray( "Depth" )->InsertTuple1( id, level );
if ( this->GenerateVectorField )
if ( this->GenerateInterfaceFields )
{
// Set vector array value
// Set interface arrays values
double v = 1. / ( 1 << level );
outData->GetArray( "Vector" )->InsertTuple3( id, v, v, v );
outData->GetArray( "Normals" )->InsertTuple3( id, v, v, v );
outData->GetArray( "Intercepts" )->InsertTuple3( id, v, 0., 3. );
}
// Initialize global index of tree
......@@ -1082,11 +1088,12 @@ void vtkHyperTreeGridSource::SubdivideFromBitsDescriptor( vtkHyperTreeGrid* outp
// Set depth array value
outData->GetArray( "Depth" )->InsertTuple1( id, level );
if ( this->GenerateVectorField )
if ( this->GenerateInterfaceFields )
{
// Set vector array value
// Set interface arrays values
double v = 1. / ( 1 << level );
outData->GetArray( "Vector" )->InsertTuple3( id, v, v, v );
outData->GetArray( "Normals" )->InsertTuple3( id, v, v, v );
outData->GetArray( "Intercepts" )->InsertTuple3( id, v, 0., 3. );
}
// Initialize global index of tree
......@@ -1303,11 +1310,12 @@ void vtkHyperTreeGridSource::SubdivideFromQuadric( vtkHyperTreeGrid* output,
// Set depth array value
outData->GetArray( "Depth" )->InsertTuple1( id, level );
if ( this->GenerateVectorField )
if ( this->GenerateInterfaceFields )
{
// Set vector array value
// Set interface arrays values
double v = 1. / ( 1 << level );
outData->GetArray( "Vector" )->InsertTuple3( id, v, v, v );
outData->GetArray( "Normals" )->InsertTuple3( id, v, v, v );
outData->GetArray( "Intercepts" )->InsertTuple3( id, v, 0., 3. );
}
// Subdivide further or stop recursion with terminal leaf
......@@ -1426,11 +1434,12 @@ void vtkHyperTreeGridSource::SubdivideFromQuadric( vtkHyperTreeGrid* output,
// Cell values
outData->GetArray( "Depth" )->InsertTuple1( id, level );
if ( this->GenerateVectorField )
if ( this->GenerateInterfaceFields )
{
// Set vector array value
// Set interface arrays values
double v = 1. / ( 1 << level );
outData->GetArray( "Vector" )->InsertTuple3( id, v, v, v );
outData->GetArray( "Normals" )->InsertTuple3( id, v, v, v );
outData->GetArray( "Intercepts" )->InsertTuple3( id, v, 0., 3. );
}
outData->GetArray( "Quadric" )->InsertTuple1( id, sum );
} // else
......
......@@ -164,13 +164,13 @@ public:
//@{
/**
* Set/get whether a cell-centered vector field
* Set/get whether cell-centered interface fields
* should be generated.
* Default: false
*/
vtkSetMacro(GenerateVectorField, bool);
vtkGetMacro(GenerateVectorField, bool);
vtkBooleanMacro(GenerateVectorField, bool);
vtkSetMacro(GenerateInterfaceFields, bool);
vtkGetMacro(GenerateInterfaceFields, bool);
vtkBooleanMacro(GenerateInterfaceFields, bool);
//@}
//@{
......@@ -327,7 +327,7 @@ protected:
unsigned int BlockSize;
bool UseDescriptor;
bool UseMaterialMask;
bool GenerateVectorField;
bool GenerateInterfaceFields;
vtkDataArray* XCoordinates;
vtkDataArray* YCoordinates;
......
Markdown is supported
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