[vtkusers] Filter to Identify Boundary Regions on any mesh shape

updega2 updega2 at berkeley.edu
Sat Jan 18 14:43:09 EST 2014


Hey Everyone,

I have been working on a filter to identify boundary regions on a mesh for
CFD purposes. I finished it up the other day, and it seems to work pretty
well. I figured I would post the code for anyone in case they happen to need
something similar in the future. Please feel free to make changes or
improvements to the code! This is my first filter, and I would love for
feedback on how to improve my code and coding practices. Thanks!

XML File:

<ServerManagerConfiguration>
    <ProxyGroup name="filters">
        <SourceProxy name="GetBoundaryFaces" class="vtkGetBoundaryFaces"
label="Boundary Face Finder">
            <Documentation
                long_help="Identify boundary surfaces with different regions
number"
                short_help="Identify boundary regions">
                The color boundary regions filter splits an external surface
in different regions. It's suitable to hidden parts of the geometry and
apply boundary conditions.
            </Documentation>
            <InputProperty
                name="Input"
                command="SetInputConnection">
                <ProxyGroupDomain name="groups">
                    <Group name="sources"/>
                    <Group name="filters"/>
                </ProxyGroupDomain>
                <DataTypeDomain name="input_type">
                    <DataType value="vtkPolyData"/>
                </DataTypeDomain>
            </InputProperty>
	    <DoubleVectorProperty command="SetFeatureAngle"
				  default_values="50.0"
				  name="FeatureAngle"
				  number_of_elements="1">
	      <DoubleRangeDomain max="180"
				 min="0"
				 name="range" />
	      <Documentation>Ths value of this property is used to define boundary
regions. If two cells share an edge boundary from a feature edges output,
then they are seperated into different regions.</Documentation>
	    </DoubleVectorProperty>
	    <Hints>
	      <Visibility replace_input="0" />
	    </Hints>
             
        </SourceProxy> 
    </ProxyGroup> 
</ServerManagerConfiguration>

vtkGetBoundaryFaces.cxx:

/*=========================================================================
 
 Program:   Visualization Toolkit
 Module:    $RCSfile: vtkFeatureEdges.cxx,v $
 
 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.
 
 =========================================================================*/
#include "vtkGetBoundaryFaces.h"

#include "vtkFloatArray.h"
#include "vtkMath.h"
#include "vtkMergePoints.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkObjectFactory.h"
#include "vtkPolyData.h"
#include "vtkPolygon.h"
#include "vtkSmartPointer.h"
#include "vtkStreamingDemandDrivenPipeline.h"
#include "vtkTriangleStrip.h"
#include "vtkUnsignedCharArray.h"
#include "vtkCellArray.h"
#include "vtkCellData.h"
#include "vtkPointData.h"
#include "vtkIncrementalPointLocator.h"
#include "vtkFeatureEdges.h"
#include "vtkCellLocator.h"

#include <iostream>

vtkCxxRevisionMacro(vtkGetBoundaryFaces, "$Revision: 0.0 $");
vtkStandardNewMacro(vtkGetBoundaryFaces);

// Construct object with feature angle = 50.0;
vtkGetBoundaryFaces::vtkGetBoundaryFaces()
{
    this->boundaries = vtkFeatureEdges::New();
    this->FeatureAngle = 50.0;
    this->newScalars = vtkIntArray::New();
    this->mesh = vtkPolyData::New();
    this->boundaryMesh = vtkPolyData::New();
    this->allIds = vtkIdList::New();
}

vtkGetBoundaryFaces::~vtkGetBoundaryFaces()
{
    if (this->boundaries)
    {
        this->boundaries->Delete();
    }
    if (this->newScalars)
    {
	this->newScalars->Delete();
    }
    if (this->mesh)
    {
	this->mesh->Delete();
    }
    if (this->boundaryMesh)
    {
	this->boundaryMesh->Delete();
    }
    if (this->allIds)
    {
	this->allIds->Delete();
    }
}

void vtkGetBoundaryFaces::PrintSelf(ostream& os, vtkIndent indent)
{
    this->Superclass::PrintSelf(os,indent);
    os << indent << "Feature Angle: " << this->FeatureAngle << "\n";
}

// Generate Separated Surfaces with Region ID Numbers
int vtkGetBoundaryFaces::RequestData(
                                 vtkInformation *vtkNotUsed(request),
                                 vtkInformationVector **inputVector,
                                 vtkInformationVector *outputVector)
{
    // get the input and output
    vtkPolyData *input = vtkPolyData::GetData(inputVector[0]);
    vtkPolyData *output = vtkPolyData::GetData(outputVector);
    
    // Define variables used by the algorithm
    int reg = 0;                                          
    vtkSmartPointer<vtkPoints> inpts = vtkSmartPointer<vtkPoints>::New();
    vtkSmartPointer<vtkCellArray> inPolys =
vtkSmartPointer<vtkCellArray>::New();
    vtkSmartPointer<vtkCellArray> inStrips =
vtkSmartPointer<vtkCellArray>::New();
    vtkSmartPointer<vtkCellArray> inLines =
vtkSmartPointer<vtkCellArray>::New();
    vtkPoints *newPts;
    vtkIdType numPts, numPolys;
    vtkIdType newId, cellId;

    //Variables for output points from Feature Edges
    vtkSmartPointer<vtkPoints> boundaryPts =
vtkSmartPointer<vtkPoints>::New();
    vtkSmartPointer<vtkCellArray> boundaryLines =
vtkSmartPointer<vtkCellArray>::New();

    //Get input points, polys and set the up in the vtkPolyData mesh
    inpts = input->GetPoints();
    inPolys = input->GetPolys();
    this->mesh->SetPoints(inpts);
    this->mesh->SetPolys(inPolys);

    //Get the number of Polys for scalar  allocation
    numPolys = input->GetNumberOfPolys();
    this->allIds->Allocate(numPolys);

    //Check the input to make sure it is there
    if (numPolys < 1)               
    {
        vtkDebugMacro("No input!");
	return 1;
    }

    //Set up Region scalar for each surface
    this->newScalars->SetNumberOfTuples(numPolys);

    //Set up Feature Edges for Boundary Edge Detection
    vtkPolyData* inputCopy = input->NewInstance();
    inputCopy->ShallowCopy(input);
    this->boundaries->SetInputData(inputCopy);
    this->boundaries->BoundaryEdgesOff();
    this->boundaries->ManifoldEdgesOff();
    this->boundaries->NonManifoldEdgesOff();
    this->boundaries->FeatureEdgesOn();
    this->boundaries->SetFeatureAngle(this->FeatureAngle);
    inputCopy->Delete();
    this->boundaries->Update();

    boundaryPts = this->boundaries->GetOutput()->GetPoints();
    boundaryLines = this->boundaries->GetOutput()->GetLines();
 
    this->boundaryMesh->SetPoints(boundaryPts);
    this->boundaryMesh->SetLines(boundaryLines);

    vtkDebugMacro("Starting Boundary Face Separation");

    //Build Links in the mesh to be able to perform complex polydata
processes;
    this->mesh->BuildLinks();

    //Set Region value of each cell to be zero initially
    for(cellId = 0; cellId < numPolys ; cellId ++)
    {
        this->newScalars->InsertValue(cellId, reg);
    }

    //Go through each cell and perfrom region identification proces
    for (cellId=0; cellId< numPolys; cellId++)
    {
       //Check to make sure the value of the region at this cellId hasn't
been set
       if (this->newScalars->GetValue(cellId) == 0)
       {
	   reg++;
	   //Call function to find all cells within certain region
	   this->FindBoundaryRegion(cellId,reg);
       }
    }

    //Copy all the input geometry and data to the output
    output->CopyStructure(input);
    output->GetPointData()->PassData(input->GetPointData());
    output->GetCellData()->PassData(input->GetCellData());

    //Add the new scalars array to the output
    this->newScalars->SetName("Regions");
    output->GetCellData()->AddArray(this->newScalars);
    output->GetCellData()->SetActiveScalars("Regions");

    return 1;
}

void vtkGetBoundaryFaces::FindBoundaryRegion(const vtkIdType cellId, int
reg)
{
    //Variables used in function
    int i;
    vtkIdType j,k,l; 
    static vtkIdType Id1=0;
    static vtkIdType Id2=0;
    vtkIdType *pts = 0;
    vtkIdType npts = 0;
    vtkIdType numNei, nei, p1, p2, nIds, neis;

    //Id List to store neighbor cells, added to member data allIds and
delted at end of function
    vtkIdList *neiIds;
    neiIds = vtkIdList::New();
    neiIds->Allocate(VTK_CELL_SIZE);
    vtkIdType numBoundPts,numBoundLines; 

    //Id List to store the one neighbor cell for each set of nodes and a
cell, deleted at end of function
    vtkIdList *neighbors;
    neighbors = vtkIdList::New();
    neighbors->Allocate(VTK_CELL_SIZE);

    //Variable for accessing neiIds list
    vtkIdType sz = 0;
    
    //Variables for the boundary cells adjacent to the boundary point
    vtkSmartPointer<vtkIdList> bLinesOne =
vtkSmartPointer<vtkIdList>::New();
    vtkSmartPointer<vtkIdList> bLinesTwo =
vtkSmartPointer<vtkIdList>::New();

    //Get output set of points and lines from the Feature Edges
    numBoundPts = this->boundaryMesh->GetPoints()->GetNumberOfPoints();
    numBoundLines = this->boundaryMesh->GetLines()->GetNumberOfCells();

    //Mark point with region number and get all points in the cell
    this->newScalars->InsertValue(cellId,reg);
    this->mesh->GetCellPoints(cellId,npts,pts);
    
    //Set up point arrays for checking of boundary points
    double boundPtComps[npts];
    double meshPt1Comps[npts];
    double meshPt2Comps[npts];

    //Get neighboring cell for each pair of points in current cell
    for (i=0; i < npts; i++)
    {
	p1 = pts[i];
	p2 = pts[(i+1)%(npts)];
	this->mesh->GetPoint(p1,meshPt1Comps);
	this->mesh->GetPoint(p2,meshPt2Comps);


	//Initial check to make sure the cell is in fact a face cell
	numNei = neighbors->GetNumberOfIds();
	this->mesh->GetCellEdgeNeighbors(cellId,p1,p2, neighbors);
	numNei = neighbors->GetNumberOfIds();

	//Check to make sure it is an oustide surface cell, i.e. one neighbor
	if (numNei==1)
	{
	    int count = 0;
	    for (j=0; j<numBoundPts;j++)
	    {

		//Check to see if cell is on the boundary, if it is get adjacent lines
		this->boundaryMesh->GetPoint(j,boundPtComps);
		if (boundPtComps[0]==meshPt1Comps[0] && boundPtComps[1]==meshPt1Comps[1]
&& boundPtComps[2]==meshPt1Comps[2])
		{ 
		    this->boundaryMesh->GetPointCells(j,bLinesOne);
		    count++;
		}
		if (boundPtComps[0]==meshPt2Comps[0] && boundPtComps[1]==meshPt2Comps[1]
&& boundPtComps[2]==meshPt2Comps[2])
		{ 
		    this->boundaryMesh->GetPointCells(j,bLinesTwo);
		    count++;
		}
	    }

	    nei=neighbors->GetId(0);
	    //if cell is not on the boundary, add new cell to check list
	    if (count < 2)
	    {
		if (this->newScalars->GetValue(nei)==0)
		{
		    neiIds->InsertId(sz,nei);
		    sz++;
		}
	    }
	    //if cell is on boundary, check to make sure it isn't false positive;
don't add to check list
	    else
	    {
		bLinesOne->IntersectWith(bLinesTwo);
		if (bLinesOne->GetNumberOfIds() == 0)
		{
		    neiIds->InsertId(sz,nei);
		    sz++;
		}
	    }
	}
    }
    
    nIds = neiIds->GetNumberOfIds();
    neighbors->Delete();
    //If there are neighboring cells in the check list, run function
recursively
    if (nIds>0)
    {
	//Add all Ids in current list to global list of Ids
	for (k=0; k< nIds;k++)
	{
	    neis = neiIds->GetId(k);
	    this->allIds->InsertId(Id1,neis);
	    Id1++;
	}
	neiIds->Delete();
	//Run function recursively for each Id just placed in the global list
	for (l=0; l< nIds;l++)
	{
	    neis = this->allIds->GetId(Id2);
            Id2++;
	    if (this->newScalars->GetValue(neis)==0)
	    {
		this->FindBoundaryRegion(neis,reg);
	    }
	}
    }
    
    else
    {
    neiIds->Delete();
    }
}

	 
vtkGetBoundaryFaces.h:

/*=========================================================================

  Program:   Visualization Toolkit
  Module:    vtkGetBoundaryFaces.h

  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.

=========================================================================*/
// .NAME vtkGetBoundaryFaces - Get Boundary Faces from poldata and label
them with integers
// .SECTION Description
// vtkGetBoundaryFaces is a filter to extract the boundary surfaces of a
model, separate the surace into multiple regions and number each region. 

// .SECTION Caveats
// To see the coloring of the liens you may have to set the ScalarMode
// instance variable of the mapper to SetScalarModeToUseCellData(). (This
// is only a problem if there are point data scalars.)

// .SECTION See Also
// vtkExtractEdges

#ifndef __vtkGetBoundaryFaces_h
#define __vtkGetBoundaryFaces_h

#include "vtkFiltersCoreModule.h" // For export macro
#include "vtkPolyDataAlgorithm.h"

class vtkFeatureEdges;

class VTKFILTERSCORE_EXPORT vtkGetBoundaryFaces : public
vtkPolyDataAlgorithm
{
public:
  static vtkGetBoundaryFaces* New();
  vtkTypeRevisionMacro(vtkGetBoundaryFaces, vtkPolyDataAlgorithm);
  void PrintSelf(ostream& os, vtkIndent indent);

  // Description:
  // Specify the feature angle for extracting feature edges.
  vtkSetClampMacro(FeatureAngle,double,0.0,180.0);
  vtkGetMacro(FeatureAngle,double);

protected:
  vtkGetBoundaryFaces();
  ~vtkGetBoundaryFaces();

  // Usual data generation method
  int RequestData(vtkInformation *vtkNotUsed(request), 
		  vtkInformationVector **inputVector, 
		  vtkInformationVector *outputVector);

  vtkFeatureEdges* boundaries;
  double FeatureAngle;
  vtkIntArray *newScalars;
  vtkPolyData *mesh;
  vtkPolyData *boundaryMesh;
  vtkIdList *allIds;

  void FindBoundaryRegion(const vtkIdType cellId, int reg);

private:
  vtkGetBoundaryFaces(const vtkGetBoundaryFaces&);  // Not implemented.
  void operator=(const vtkGetBoundaryFaces&);  // Not implemented.
};

#endif






--
View this message in context: http://vtk.1045678.n5.nabble.com/Filter-to-Identify-Boundary-Regions-on-any-mesh-shape-tp5725505.html
Sent from the VTK - Users mailing list archive at Nabble.com.


More information about the vtkusers mailing list