[vtkusers] reading element blocks from exodus mesh into a vtkSurfaceFilter

[Joseph Butner]

Hi all,


I'm going to answer my own question in the hopes that it helps someone else.  Note that the exodus ii format does not allow element block indices of 0 (at least when designing meshes with Trelis/Cubit), so if your mesh has indices  [1, 2....n-1, n]  they will correspond to indices [0, 1........n-2, n-1] in the std::vector.


The following will read numbered element block node sets (in numerical order) and pass them to a vtkSurfaceFilter (in this case, a vector of surface filters):


#define Instantiate(obj, class) vtkSmartPointer<class> obj = vtkSmartPointer<class>::New()


    Instantiate( exodusiiReader, vtkExodusIIReader );
    std::vector<vtkSmartPointer<vtkPolyData> > element_blocks_surface_filter_vector;


    std::string exodusMesh = "/path/to/mesh/mesh_name.exo";
    exodusiiReader->SetFileName(exodusMesh.c_str());
    exodusiiReader->Update();
    exodusiiReader->UpdateInformation();  // read mesh metadata

    vtkCompositeDataIterator* iter = (exodusiiReader->GetOutput())->NewIterator();
    iter->GoToFirstItem();
    while (!iter->IsDoneWithTraversal())
    {
          vtkDataObject* dObj = iter->GetCurrentDataObject();

            Instantiate( surfaceFilter,vtkDataSetSurfaceFilter );
        #if VTK_MAJOR_VERSION <= 5
          surfaceFilter->SetInput(dObj);
        #else
          surfaceFilter->SetInputData(dObj);
        #endif
          surfaceFilter->Update();

          element_blocks_surface_filter_vector.push_back(surfaceFilter->GetOutput());

          iter->GoToNextItem();
    }


I hope this makes it to the net and helps others in the future.


Joe


________________________________
From: vtkusers <vtkusers-bounces at vtk.org> on behalf of Joseph Butner <liljoe at unm.edu>
Sent: Friday, December 12, 2014 11:47 AM
To: vtkusers at vtk.org
Subject: [vtkusers] reading element blocks from exodus mesh into a vtkSurfaceFilter


I have an exodus mesh file that contains numbered element blocks.  I want to be able to determine which element block (if any) any (x, y, z) coordinate is in.  To do this, I read the mesh in with vtkExodusiiMeshReader, attempt to generate a vtkSurfaceFilter for each element block, and use vtkSelectEnclosedPoints to find the element block the point is in.  Unfortunately, I'm doing something wrong when reading the nodes for each element block, the data structure that is supposed to have the node coordinates seems to be empty, it throws this error:

ERROR: In /home/joe/Downloads/VTK-6.1.0/Common/ExecutionModel/vtkDemandDrivenPipeline.cxx, line 715
vtkCompositeDataPipeline (0x35181a0): Input port 1 of algorithm vtkSelectEnclosedPoints(0x1c31480) has 0 connections but is not optional.


A small piece of the code I have is:


#define Instantiate(obj, class) vtkSmartPointer<class> obj = vtkSmartPointer<class>::New()

//I want to store the vtkPolyData for each element block in a vector, will be passed to vtkSelectEnclosedPoints later:
std::vector<vtkSmartPointer<vtkPolyData> > element_blocks_surface_filter_vector;


std::string exodusMesh = "foo.exo";
exodusiiReader->SetFileName(exodusMesh.c_str());
exodusiiReader->Update();


// get the metadata from the mesh:
exodusiiReader->UpdateInformation();

// read in the element blocks, make a surface filter for each, push_back into vector:
for (int z = 0; z < exodusiiReader->GetNumberOfElementBlockArrays(); z++)
{
    const char* nname = exodusiiReader->GetElementBlockArrayName(z);

/*
 *  I think my mistake is in the following line
 */
    exodusiiReader->SetObjectArrayStatus(exodusiiReader->ELEM_BLOCK, nname, 1);
    exodusiiReader->Update();

    Instantiate( surfaceFilter,vtkDataSetSurfaceFilter );
    Instantiate( MultiBlockData, vtkMultiBlockDataSet );
    Instantiate( DataObject, vtkDataObject );

    MultiBlockData = exodusiiReader->GetOutput();
    DataObject= MultiBlockData->GetBlock(0);


/*
 * everything below here should be fine (except the DataObject generated above isn't correct).

 * it works when I just send it the whole mesh
 */

    #if VTK_MAJOR_VERSION <= 5
      surfaceFilter->SetInput(DataObject);
    #else
      surfaceFilter->SetInputData(DataObject);
    #endif
      surfaceFilter->Update();

    element_blocks_surface_filter_vector.push_back(surfaceFilter->GetOutput());

}


Any information about what I'm doing wrong above would be greatly appreciated. Thanks!


Joe
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/vtkusers/attachments/20141217/3b66f4fa/attachment.html>