vtkSubdivideTetra.cxx 6.61 KB
Newer Older
1 2 3 4 5
/*=========================================================================

  Program:   Visualization Toolkit
  Module:    vtkSubdivideTetra.cxx

6
  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
7 8
  All rights reserved.
  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
9

10 11
     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
12
     PURPOSE.  See the above copyright notice for more information.
13 14 15

=========================================================================*/
#include "vtkSubdivideTetra.h"
16

17 18
#include "vtkCellType.h"
#include "vtkGenericCell.h"
19 20
#include "vtkInformation.h"
#include "vtkInformationVector.h"
21
#include "vtkMergePoints.h"
22
#include "vtkObjectFactory.h"
23
#include "vtkPointData.h"
24
#include "vtkUnstructuredGrid.h"
25

Brad King's avatar
Brad King committed
26
vtkStandardNewMacro(vtkSubdivideTetra);
27

28
//----------------------------------------------------------------------------
29 30 31 32 33 34
// Description:
// Construct with all types of clipping turned off.
vtkSubdivideTetra::vtkSubdivideTetra()
{
}

35
//----------------------------------------------------------------------------
36 37 38 39
int vtkSubdivideTetra::RequestData(
  vtkInformation *vtkNotUsed(request),
  vtkInformationVector **inputVector,
  vtkInformationVector *outputVector)
40
{
41 42 43 44
  // get the info objects
  vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
  vtkInformation *outInfo = outputVector->GetInformationObject(0);

45
  // get the input and output
46 47 48 49 50
  vtkUnstructuredGrid *input = vtkUnstructuredGrid::SafeDownCast(
    inInfo->Get(vtkDataObject::DATA_OBJECT()));
  vtkUnstructuredGrid *output = vtkUnstructuredGrid::SafeDownCast(
    outInfo->Get(vtkDataObject::DATA_OBJECT()));

Amy Squillacote's avatar
Amy Squillacote committed
51 52
  vtkIdType numPts = input->GetNumberOfPoints();
  vtkIdType numCells = input->GetNumberOfCells();
53
  vtkPoints *inPts = input->GetPoints();
Amy Squillacote's avatar
Amy Squillacote committed
54
  vtkIdType cellId, i;
55
  vtkIdType pts[4];
56
  vtkGenericCell *cell;
57 58 59
  vtkPointData *pd = input->GetPointData();
  vtkPointData *outputPD = output->GetPointData();
  vtkPoints *newPts;
Amy Squillacote's avatar
Amy Squillacote committed
60
  vtkIdType ptId;
61

Ken Martin's avatar
Ken Martin committed
62
  double weights[4], x0[3], x1[3], x2[3], x3[3], x[3];
Amy Squillacote's avatar
Amy Squillacote committed
63 64
  int p0, p1, p2, p3;
  vtkIdType center, e01, e02, e03, e12, e13, e23;
65
  vtkMergePoints *locator;
66

67 68
  vtkDebugMacro(<<"Executing mesh subdivide");

Jeff Lee's avatar
Jeff Lee committed
69
  if (input->IsHomogeneous() == 0 ||
70
      input->GetCellType(0) != VTK_TETRA)
71
    {
72
    vtkErrorMacro(<<"all cells must be tetrahedra.");
73
    return 1;
74
    }
75

76 77 78 79
  // Copy original points and point data
  newPts = vtkPoints::New();
  newPts->Allocate(5*numPts,numPts);
  outputPD->InterpolateAllocate(pd,5*numPts,numPts);
80

81 82 83 84 85 86 87 88 89 90 91 92
  output->Allocate(numCells);
  output->SetPoints(newPts);

  locator = vtkMergePoints::New();
  locator->InitPointInsertion (newPts, input->GetBounds());

  for (ptId=0; ptId < numPts; ptId++)
    {
    locator->InsertNextPoint(inPts->GetPoint(ptId));
    outputPD->CopyData(pd,ptId,ptId);
    }

93 94
  cell = vtkGenericCell::New();

95 96 97 98
  // loop over tetrahedra, generating sixteen new ones for each. This is
  // done by introducing mid-edge nodes and a single mid-tetra node.
  for(cellId=0; cellId < numCells; cellId++)
    {
99
    input->GetCell(cellId, cell);
100 101

    // get tetra points
102 103 104 105
    cell->Points->GetPoint(0,x0);
    cell->Points->GetPoint(1,x1);
    cell->Points->GetPoint(2,x2);
    cell->Points->GetPoint(3,x3);
106

107 108 109 110
    p0 = cell->PointIds->GetId(0);
    p1 = cell->PointIds->GetId(1);
    p2 = cell->PointIds->GetId(2);
    p3 = cell->PointIds->GetId(3);
111

112
    // compute center point
113
    weights[0] = weights[1] = weights[2] = weights[3] = 0.25;
Bill Lorensen's avatar
Bill Lorensen committed
114 115 116 117
    for (i=0; i<3; i++)
      {
      x[i] = 0.25*(x0[i] + x1[i] + x2[i] + x3[i]);
      }
118
    center = locator->InsertNextPoint(x);
119
    outputPD->InterpolatePoint(pd, center, cell->PointIds, weights);
120

121 122
    // compute edge points
    // edge 0-1
Bill Lorensen's avatar
Bill Lorensen committed
123 124 125 126
    for (i=0; i<3; i++)
      {
      x[i] = 0.5 * (x1[i] + x0[i]);
      }
127 128
    e01 = locator->InsertNextPoint(x);
    outputPD->InterpolateEdge(pd, e01, p0, p1, 0.5);
129

130
    // edge 1-2
Bill Lorensen's avatar
Bill Lorensen committed
131 132 133 134
    for (i=0; i<3; i++)
      {
      x[i] = 0.5 * (x2[i] + x1[i]);
      }
135 136
    e12 = locator->InsertNextPoint(x);
    outputPD->InterpolateEdge(pd, e12, p1, p2, 0.5);
137

138
    // edge 2-0
Bill Lorensen's avatar
Bill Lorensen committed
139 140 141 142
    for (i=0; i<3; i++)
      {
      x[i] = 0.5 * (x2[i] + x0[i]);
      }
143 144
    e02 = locator->InsertNextPoint(x);
    outputPD->InterpolateEdge(pd, e02, p2, p0, 0.5);
145

146
    // edge 0-3
Bill Lorensen's avatar
Bill Lorensen committed
147 148 149 150
    for (i=0; i<3; i++)
      {
      x[i] = 0.5 * (x3[i] + x0[i]);
      }
151 152
    e03 = locator->InsertNextPoint(x);
    outputPD->InterpolateEdge(pd, e03, p0, p3, 0.5);
153

154
    // edge 1-3
Bill Lorensen's avatar
Bill Lorensen committed
155 156 157 158
    for (i=0; i<3; i++)
      {
      x[i] = 0.5 * (x3[i] + x1[i]);
      }
159 160
    e13 = locator->InsertNextPoint(x);
    outputPD->InterpolateEdge(pd, e13, p1, p3, 0.5);
161

162
    // edge 2-3
Bill Lorensen's avatar
Bill Lorensen committed
163 164 165 166
    for (i=0; i<3; i++)
      {
      x[i] = 0.5 * (x3[i] + x2[i]);
      }
167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231
    e23 = locator->InsertNextPoint(x);
    outputPD->InterpolateEdge(pd, e23, p2, p3, 0.5);

    // Now create tetrahedra
    // First, four tetra from each vertex
    pts[0] = p0;
    pts[1] = e01;
    pts[2] = e02;
    pts[3] = e03;
    output->InsertNextCell(VTK_TETRA, 4, pts);
    pts[0] = p1;
    pts[1] = e01;
    pts[2] = e12;
    pts[3] = e13;
    output->InsertNextCell(VTK_TETRA, 4, pts);
    pts[0] = p2;
    pts[1] = e02;
    pts[2] = e12;
    pts[3] = e23;
    output->InsertNextCell(VTK_TETRA, 4, pts);
    pts[0] = p3;
    pts[1] = e03;
    pts[2] = e13;
    pts[3] = e23;
    output->InsertNextCell(VTK_TETRA, 4, pts);

    // Now four tetra from cut-off tetra corners
    pts[0] = center;
    pts[1] = e01;
    pts[2] = e02;
    pts[3] = e03;
    output->InsertNextCell(VTK_TETRA, 4, pts);
    pts[1] = e01;
    pts[2] = e12;
    pts[3] = e13;
    output->InsertNextCell(VTK_TETRA, 4, pts);
    pts[1] = e02;
    pts[2] = e12;
    pts[3] = e23;
    output->InsertNextCell(VTK_TETRA, 4, pts);
    pts[1] = e03;
    pts[2] = e13;
    pts[3] = e23;
    output->InsertNextCell(VTK_TETRA, 4, pts);

    // Now four tetra from triangles on tetra faces
    pts[0] = center;
    pts[1] = e01;
    pts[2] = e12;
    pts[3] = e02;
    output->InsertNextCell(VTK_TETRA, 4, pts);
    pts[1] = e01;
    pts[2] = e13;
    pts[3] = e03;
    output->InsertNextCell(VTK_TETRA, 4, pts);
    pts[1] = e12;
    pts[2] = e23;
    pts[3] = e13;
    output->InsertNextCell(VTK_TETRA, 4, pts);
    pts[1] = e02;
    pts[2] = e23;
    pts[3] = e03;
    output->InsertNextCell(VTK_TETRA, 4, pts);

    } //for all cells
Will Schroeder's avatar
Will Schroeder committed
232
  cell->Delete();
233 234 235

  vtkDebugMacro(<<"Subdivided " << numCells << " cells");

236
  locator->Delete();
237 238
  newPts->Delete();
  output->Squeeze();
239 240

  return 1;
241 242
}

243
//----------------------------------------------------------------------------
244
void vtkSubdivideTetra::PrintSelf(ostream& os, vtkIndent indent)
245
{
Brad King's avatar
Brad King committed
246
  this->Superclass::PrintSelf(os,indent);
247
}