Merge topic 'AddYoungsMaterialInterfaceReconstuctionExample'

8fdead8b A much better example using blocs for per-material meshes
a98720c4 More generically add interface mappers/actors with composite iterator
c7b0951c Improvements to the YMI class
291c30d8 Show reconstructed interface
d08cd68c Added Youngs material interface instance along with mapper/actor
ca6d081f Wireframe overlay of presence fraction #1
4a878122 An example reading AVS/UCD binary file for Youngs Material Interface

Examples/Graphics/Python/YoungsMaterialInterface.py

+############################################################
+from vtk import *
+from vtk.util.misc import vtkGetDataRoot
+############################################################
+
+# Create renderer and add actors to it
+renderer = vtkRenderer()
+renderer.SetBackground( .8, .8 ,.8 )
+
+# Create render window
+window = vtkRenderWindow()
+window.AddRenderer( renderer )
+window.SetSize( 500, 500 )
+
+# Create interactor
+interactor = vtkRenderWindowInteractor()
+interactor.SetRenderWindow( window )
+
+# Retrieve data root
+VTK_DATA_ROOT = vtkGetDataRoot()
+
+# Read from AVS UCD data in binary form
+reader = vtkAVSucdReader()
+reader.SetFileName( VTK_DATA_ROOT + "/Data/UCD/UCD_00005.inp" )
+
+# Update reader and get mesh cell data
+reader.Update()
+mesh = reader.GetOutput()
+cellData = mesh.GetCellData()
+
+# Create normal vectors
+cellData.SetActiveScalars( "norme[0]" )
+normX = cellData.GetScalars()
+cellData.SetActiveScalars( "norme[1]" )
+normY = cellData.GetScalars()
+n = normX.GetNumberOfTuples()
+norm = vtkDoubleArray()
+norm.SetNumberOfComponents( 3 )
+norm.SetNumberOfTuples( n )
+norm.SetName( "norme" )
+for i in range( 0, n ):
+    norm.SetTuple3( i, normX.GetTuple1( i ), normY.GetTuple1( i ), 0.0 )
+cellData.SetVectors( norm )
+
+# Extract submeshes corresponding to 2 different material Ids
+cellData.SetActiveScalars( "Material Id" )
+threshold2 = vtkThreshold()
+threshold2.SetInput( mesh )
+threshold2.SetInputArrayToProcess(0, 0, 0, vtkDataObject.FIELD_ASSOCIATION_CELLS, vtkDataSetAttributes.SCALARS );
+threshold2.ThresholdByLower( 2 )
+threshold2.Update()
+meshMat2 = threshold2.GetOutput()
+threshold3 = vtkThreshold()
+threshold3.SetInput( mesh )
+threshold3.SetInputArrayToProcess(0, 0, 0, vtkDataObject.FIELD_ASSOCIATION_CELLS, vtkDataSetAttributes.SCALARS );
+threshold3.ThresholdByUpper( 3 )
+threshold3.Update()
+meshMat3 = threshold3.GetOutput()
+
+# Make multiblock from extracted submeshes
+meshMB =  vtkMultiBlockDataSet()
+meshMB.SetNumberOfBlocks( 2 )
+meshMB.GetMetaData( 0 ).Set( vtkCompositeDataSet.NAME(), "Material 2" )
+meshMB.SetBlock( 0, meshMat2 )
+meshMB.GetMetaData( 1 ).Set( vtkCompositeDataSet.NAME(), "Material 3" )
+meshMB.SetBlock( 1, meshMat3 )
+
+# Create mapper for both submeshes
+matRange = cellData.GetScalars().GetRange()
+mat2Mapper = vtkDataSetMapper()
+mat2Mapper.SetInput( meshMat2 )
+mat2Mapper.SetScalarRange( matRange )
+mat2Mapper.SetScalarModeToUseCellData()
+mat2Mapper.SetColorModeToMapScalars()
+mat2Mapper.ScalarVisibilityOn()
+mat2Mapper.SetResolveCoincidentTopologyPolygonOffsetParameters( 0, 1 )
+mat2Mapper.SetResolveCoincidentTopologyToPolygonOffset()
+mat3Mapper = vtkDataSetMapper()
+mat3Mapper.SetInput( meshMat3 )
+mat3Mapper.SetScalarRange( matRange )
+mat3Mapper.SetScalarModeToUseCellData()
+mat3Mapper.SetColorModeToMapScalars()
+mat3Mapper.ScalarVisibilityOn()
+mat3Mapper.SetResolveCoincidentTopologyPolygonOffsetParameters( 1, 1 )
+mat3Mapper.SetResolveCoincidentTopologyToPolygonOffset()
+
+# Create wireframe actor for one, surface actor for other, and add both to view
+wireActor = vtkActor()
+wireActor.SetMapper( mat2Mapper )
+wireActor.GetProperty().SetRepresentationToWireframe()
+renderer.AddViewProp( wireActor )
+surfActor = vtkActor()
+surfActor.SetMapper( mat3Mapper )
+surfActor.GetProperty().SetRepresentationToSurface()
+renderer.AddViewProp( surfActor )
+
+cellData.SetActiveScalars("frac_pres[1]")
+# Reconstruct material interface
+interface = vtkYoungsMaterialInterface()
+interface.SetInput( meshMB )
+interface.SetNumberOfMaterials( 2 )
+interface.SetMaterialVolumeFractionArray( 0, "frac_pres[1]" )
+interface.SetMaterialVolumeFractionArray( 1, "frac_pres[2]" )
+interface.SetMaterialNormalArray( 0, "norme" )
+interface.SetMaterialNormalArray( 1, "norme" )
+interface.UseAllBlocksOn()
+interface.Update()
+
+# Create mappers and actors for surface rendering of all reconstructed interfaces
+interfaceIterator = vtkCompositeDataIterator()
+interfaceIterator.SetDataSet( interface.GetOutput() )
+interfaceIterator.VisitOnlyLeavesOn()
+interfaceIterator.SkipEmptyNodesOn()
+interfaceIterator.InitTraversal()
+interfaceIterator.GoToFirstItem()
+while ( interfaceIterator.IsDoneWithTraversal() == 0 ):
+    idx = interfaceIterator.GetCurrentFlatIndex()
+    # Create mapper for leaf node
+    print "Creating mapper and actor for object with flat index", idx
+    interfaceMapper = vtkDataSetMapper()
+    interfaceMapper.SetInput( interfaceIterator.GetCurrentDataObject() )
+    interfaceIterator.GoToNextItem()
+    interfaceMapper.ScalarVisibilityOff()
+    # Create surface actor and add it to view
+    interfaceActor = vtkActor()
+    interfaceActor.SetMapper( interfaceMapper )
+    if ( idx == 2 ):
+        interfaceActor.GetProperty().SetColor( 0, 1, 0 )
+    else:
+        interfaceActor.GetProperty().SetColor( 0, 0, 0 )
+    renderer.AddViewProp( interfaceActor )
+
+# Start interaction
+window.Render()
+interactor.Start()

Graphics/vtkYoungsMaterialInterface.cxx

@@ -453,7 +453,7 @@ void vtkYoungsMaterialInterface::UpdateBlockMapping()
 {
   int n = this->MaterialBlockMapping->GetNumberOfTuples();
   int curmat = -1;
-  for(int i = 0;i<n;i++)
+  for( int i = 0; i < n; ++ i )
     {
     int b = this->MaterialBlockMapping->GetValue(i);
     vtkDebugMacro(<<"MaterialBlockMapping "<<b<<"\n");
@@ -517,7 +517,8 @@ int vtkYoungsMaterialInterface::RequestData(
   vtkIdType debugStats_NullNormal = 0;
   vtkIdType debugStats_NoInterfaceFound = 0;
 
-  int nmat = static_cast<int> (this->Internals->Materials.size());
+  // Initialize number of materials
+  int nmat = static_cast<int>( this->Internals->Materials.size() );
 
   // alocate composite iterator
   vtkSmartPointer<vtkCompositeDataIterator> inputIterator = vtkSmartPointer<vtkCompositeDataIterator>::New();
@@ -529,17 +530,20 @@ int vtkYoungsMaterialInterface::RequestData(
 
   // first compute number of domains
   int* inputsPerMaterial = new int[nmat];
-  for (int i = 0; i < nmat; i++) { inputsPerMaterial[i] = 0; }
+  for ( int i = 0; i < nmat; ++ i ) 
+    {
+    inputsPerMaterial[i] = 0;
+    }
 
-  while (inputIterator->IsDoneWithTraversal() == 0)
+  while ( ! inputIterator->IsDoneWithTraversal() )
     {
-    vtkDataSet * input = vtkDataSet::SafeDownCast(
-                                                  inputIterator->GetCurrentDataObject());
+    vtkDataSet * input 
+      = vtkDataSet::SafeDownCast( inputIterator->GetCurrentDataObject() );
     // Composite indices begin at 1 (0 is the root)
     int composite_index = inputIterator->GetCurrentFlatIndex();
     inputIterator->GoToNextItem();
 
-    if (input != 0 && input->GetNumberOfCells() > 0)
+    if ( input && input->GetNumberOfCells() > 0 )
       {
       int m = 0;
       for (std::vector<vtkYoungsMaterialInterfaceInternals::MaterialDescription>::iterator
@@ -547,15 +551,15 @@ int vtkYoungsMaterialInterface::RequestData(
            it!= this->Internals->Materials.end(); ++it, ++m)
         {
         vtkDataArray* fraction = input->GetCellData()->GetArray((*it).volume.c_str());
-        bool materialHasBlock = ((*it).blocks.find(composite_index)!= (*it).blocks.end());
+        bool materialHasBlock = ( (*it).blocks.find(composite_index)!= (*it).blocks.end() );
         if ( fraction && 
-             ( this->UseAllBlocks ||materialHasBlock ) )
+             ( this->UseAllBlocks || materialHasBlock ) )
           {
           double range[2];
           fraction->GetRange(range);
           if (range[1] > this->VolumeFractionRange[0])
             {
-            ++inputsPerMaterial[m];
+            ++ inputsPerMaterial[m];
             }
           }
         }
@@ -616,35 +620,39 @@ int vtkYoungsMaterialInterface::RequestData(
              it = this->Internals->Materials.begin(); it
              != this->Internals->Materials.end(); ++it, ++m)
         {
-        Mats[m].fractionArray = input->GetCellData()->GetArray(
-                                                               (*it).volume.c_str());
-        Mats[m].normalArray = input->GetCellData()->GetArray(
-                                                             (*it).normal.c_str());
-        Mats[m].normalXArray = input->GetCellData()->GetArray(
-                                                              (*it).normalX.c_str());
-        Mats[m].normalYArray = input->GetCellData()->GetArray(
-                                                              (*it).normalY.c_str());
-        Mats[m].normalZArray = input->GetCellData()->GetArray(
-                                                              (*it).normalZ.c_str());
-        Mats[m].orderingArray = input->GetCellData()->GetArray(
-                                                               (*it).ordering.c_str());
-
-        if (Mats[m].fractionArray == 0)
+        Mats[m].fractionArray 
+          = input->GetCellData()->GetArray( (*it).volume.c_str() );
+        Mats[m].normalArray 
+          = input->GetCellData()->GetArray( (*it).normal.c_str() );
+        Mats[m].normalXArray 
+          = input->GetCellData()->GetArray( (*it).normalX.c_str() );
+        Mats[m].normalYArray 
+          = input->GetCellData()->GetArray( (*it).normalY.c_str() );
+        Mats[m].normalZArray 
+          = input->GetCellData()->GetArray( (*it).normalZ.c_str() );
+        Mats[m].orderingArray 
+          = input->GetCellData()->GetArray( (*it).ordering.c_str() );
+
+        if ( ! Mats[m].fractionArray )
           {
           vtkDebugMacro(<<"Material "<<m<<": volume fraction array '"<<(*it).volume<<"' not found\n");
           }
-        //                      if( Mats[m].orderingArray == 0 )
-        //                      {
-        //                              vtkDebugMacro(<<"Material "<<m<<" material ordering array '"<<(*it).ordering<<"' not found\n");
-        //                      }
-        if( Mats[m].normalArray==0 && Mats[m].normalXArray==0 && Mats[m].normalYArray==0 && Mats[m].normalZArray==0 )
+        if( ! Mats[m].orderingArray )
+          {
+          vtkDebugMacro(<<"Material "<<m<<" material ordering array '"<<(*it).ordering<<"' not found\n");
+          }
+        if( ! Mats[m].normalArray 
+            && ! Mats[m].normalXArray 
+            && ! Mats[m].normalYArray 
+            && ! Mats[m].normalZArray )
           {
           vtkDebugMacro(<<"Material "<<m<<" normal  array '"<<(*it).normal<<"' not found\n");
           }
 
         bool materialHasBlock = ( (*it).blocks.find(composite_index) != (*it).blocks.end() );
-        if( !materialHasBlock )
+        if( ! this->UseAllBlocks && ! materialHasBlock )
           {
+          // FIXME: verify the former condition
           Mats[m].fractionArray = 0; // TODO: we certainly can do better to avoid material calculations
           }
 
@@ -653,7 +661,7 @@ int vtkYoungsMaterialInterface::RequestData(
         Mats[m].cellArrayCount = 0;
 
         Mats[m].outCellArrays = new vtkDataArray* [ nCellData ];
-        for(int i = 0;i<nCellData;i++)
+        for( int i = 0; i < nCellData; ++ i )
           {
           Mats[m].outCellArrays[i] = vtkDataArray::CreateDataArray( inCellArrays[i]->GetDataType() );
           Mats[m].outCellArrays[i]->SetName( inCellArrays[i]->GetName() );
@@ -664,7 +672,7 @@ int vtkYoungsMaterialInterface::RequestData(
         Mats[m].pointCount = 0;
         Mats[m].outPointArrays = new vtkDataArray* [ nPointData ];
 
-        for(int i = 0;i<(nPointData-1);i++)
+        for( int i = 0;i<(nPointData-1);i++)
           {
           Mats[m].outPointArrays[i] = vtkDataArray::CreateDataArray( inPointArrays[i]->GetDataType() );
           Mats[m].outPointArrays[i]->SetName( inPointArrays[i]->GetName() );