1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkDataSetGradient.cxx
5
6 Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
7 All rights reserved.
8 See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
9
10 This software is distributed WITHOUT ANY WARRANTY; without even
11 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
12 PURPOSE. See the above copyright notice for more information.
13
14 =========================================================================*/
15 // .SECTION Thanks
16 // This file is part of the generalized Youngs material interface reconstruction algorithm
17 // contributed by CEA/DIF - Commissariat a l'Energie Atomique, Centre DAM Ile-De-France <br> BP12,
18 // F-91297 Arpajon, France. <br> Implementation by Thierry Carrard (CEA)
19
20 #include "vtkDataSetGradient.h"
21 #include "vtkCell.h"
22 #include "vtkCellData.h"
23 #include "vtkDataArray.h"
24 #include "vtkDataSet.h"
25 #include "vtkDataSetGradientPrecompute.h"
26 #include "vtkDoubleArray.h"
27 #include "vtkInformation.h"
28 #include "vtkInformationVector.h"
29 #include "vtkMath.h"
30 #include "vtkObjectFactory.h"
31 #include "vtkPointData.h"
32
33 // utility macros
34 #define ADD_VEC(a, b) \
35 a[0] += b[0]; \
36 a[1] += b[1]; \
37 a[2] += b[2]
38 #define SCALE_VEC(a, b) \
39 a[0] *= b; \
40 a[1] *= b; \
41 a[2] *= b
42 #define ZERO_VEC(a) \
43 a[0] = 0; \
44 a[1] = 0; \
45 a[2] = 0
46 #define COPY_VEC(a, b) \
47 a[0] = b[0]; \
48 a[1] = b[1]; \
49 a[2] = b[2];
50 #define MAX_CELL_POINTS 128
51 #define VTK_CQS_EPSILON 1e-12
52
53 // standard constructors and factory
54 vtkStandardNewMacro(vtkDataSetGradient);
55
56 /*!
57 The default constructor
58 \sa ~vtkDataSetGradient()
59 */
vtkDataSetGradient()60 vtkDataSetGradient::vtkDataSetGradient()
61 : ResultArrayName(nullptr)
62 {
63 this->SetResultArrayName("gradient");
64 }
65
66 /*!
67 The destructor
68 \sa vtkDataSetGradient()
69 */
~vtkDataSetGradient()70 vtkDataSetGradient::~vtkDataSetGradient()
71 {
72 this->SetResultArrayName(nullptr);
73 }
74
PrintSelf(ostream & os,vtkIndent indent)75 void vtkDataSetGradient::PrintSelf(ostream& os, vtkIndent indent)
76 {
77 this->Superclass::PrintSelf(os, indent);
78 os << indent << "Result array name: " << this->ResultArrayName << "\n";
79 }
80
RequestData(vtkInformation * vtkNotUsed (request),vtkInformationVector ** inputVector,vtkInformationVector * outputVector)81 int vtkDataSetGradient::RequestData(vtkInformation* vtkNotUsed(request),
82 vtkInformationVector** inputVector, vtkInformationVector* outputVector)
83 {
84 // get the info objects
85 vtkInformation* inInfo = inputVector[0]->GetInformationObject(0);
86 vtkInformation* outInfo = outputVector->GetInformationObject(0);
87
88 // get connected input & output
89 vtkDataSet* _output = vtkDataSet::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT()));
90 vtkDataSet* _input = vtkDataSet::SafeDownCast(inInfo->Get(vtkDataObject::DATA_OBJECT()));
91
92 if (_input == nullptr || _output == nullptr)
93 {
94 vtkErrorMacro(<< "Missing input or output \n");
95 return 0;
96 }
97
98 // get array to compute gradient from
99 vtkDataArray* inArray = this->GetInputArrayToProcess(0, _input);
100 if (inArray == nullptr)
101 {
102 inArray = _input->GetPointData()->GetScalars();
103 }
104 if (inArray == nullptr)
105 {
106 inArray = _input->GetCellData()->GetScalars();
107 }
108
109 if (inArray == nullptr)
110 {
111 vtkErrorMacro(<< "no input array to process\n");
112 return 0;
113 }
114
115 vtkDebugMacro(<< "Input array to process : " << inArray->GetName() << "\n");
116
117 bool pointData;
118 if (_input->GetCellData()->GetArray(inArray->GetName()) == inArray)
119 {
120 pointData = false;
121 vtkDebugMacro(<< "cell data to point gradient\n");
122 }
123 else if (_input->GetPointData()->GetArray(inArray->GetName()) == inArray)
124 {
125 pointData = true;
126 vtkDebugMacro(<< "point data to cell gradient\n");
127 }
128 else
129 {
130 vtkErrorMacro(<< "input array must be cell or point data\n");
131 return 0;
132 }
133
134 // we're just adding a scalar field
135 _output->ShallowCopy(_input);
136
137 vtkDataArray* cqsArray = _output->GetFieldData()->GetArray("GradientPrecomputation");
138 vtkDataArray* sizeArray = _output->GetCellData()->GetArray("CellSize");
139 if (cqsArray == nullptr || sizeArray == nullptr)
140 {
141 vtkDebugMacro(
142 << "Couldn't find field array 'GradientPrecomputation', computing it right now.\n");
143 vtkDataSetGradientPrecompute::GradientPrecompute(_output);
144 cqsArray = _output->GetFieldData()->GetArray("GradientPrecomputation");
145 sizeArray = _output->GetCellData()->GetArray("CellSize");
146 if (cqsArray == nullptr || sizeArray == nullptr)
147 {
148 vtkErrorMacro(
149 << "Computation of field array 'GradientPrecomputation' or 'CellSize' failed.\n");
150 return 0;
151 }
152 }
153
154 vtkIdType nCells = _input->GetNumberOfCells();
155 vtkIdType nPoints = _input->GetNumberOfPoints();
156
157 vtkDoubleArray* gradientArray = vtkDoubleArray::New();
158 gradientArray->SetName(this->ResultArrayName);
159 gradientArray->SetNumberOfComponents(3);
160
161 if (pointData) // compute cell gradient from point data
162 {
163 gradientArray->SetNumberOfTuples(nCells);
164 vtkIdType cellPoint = 0;
165 for (vtkIdType i = 0; i < nCells; i++)
166 {
167 vtkCell* cell = _input->GetCell(i);
168 int np = cell->GetNumberOfPoints();
169 double gradient[3] = { 0, 0, 0 };
170 for (int p = 0; p < np; p++)
171 {
172 double cqs[3], scalar;
173 cqsArray->GetTuple(cellPoint++, cqs);
174 scalar = inArray->GetTuple1(cell->GetPointId(p));
175 SCALE_VEC(cqs, scalar);
176 ADD_VEC(gradient, cqs);
177 }
178 SCALE_VEC(gradient, (1.0 / sizeArray->GetTuple1(i)));
179 gradientArray->SetTuple(i, gradient);
180 }
181
182 _output->GetCellData()->AddArray(gradientArray);
183 //_output->GetCellData()->SetVectors( gradientArray );
184 }
185 else // compute point gradient from cell data
186 {
187 gradientArray->SetNumberOfTuples(nPoints);
188 gradientArray->FillComponent(0, 0.0);
189 gradientArray->FillComponent(1, 0.0);
190 gradientArray->FillComponent(2, 0.0);
191 double* gradient = gradientArray->WritePointer(0, nPoints * 3);
192 double* gradientDivisor = new double[nPoints];
193 for (vtkIdType i = 0; i < nPoints; i++)
194 {
195 gradientDivisor[i] = 0.0;
196 }
197 vtkIdType cellPoint = 0;
198 for (vtkIdType i = 0; i < nCells; i++)
199 {
200 vtkCell* cell = _input->GetCell(i);
201 int np = cell->GetNumberOfPoints();
202 double scalar = inArray->GetTuple1(i);
203 for (int p = 0; p < np; p++)
204 {
205 double cqs[3];
206 double pointCoord[3];
207 vtkIdType pointId = cell->GetPointId(p);
208 cqsArray->GetTuple(cellPoint++, cqs);
209 _input->GetPoint(cell->GetPointId(p), pointCoord);
210 scalar *= cell->GetCellDimension();
211 SCALE_VEC((gradient + pointId * 3), scalar);
212 gradientDivisor[pointId] += vtkMath::Dot(cqs, pointCoord);
213 }
214 }
215 for (vtkIdType i = 0; i < nPoints; i++)
216 {
217 SCALE_VEC((gradient + i * 3), (1.0 / gradientDivisor[i]));
218 }
219 delete[] gradientDivisor;
220 _output->GetPointData()->AddArray(gradientArray);
221 //_output->GetPointData()->SetVectors( gradientArray );
222 }
223 gradientArray->Delete();
224
225 vtkDebugMacro(<< _output->GetClassName() << " @ " << _output << " :\n");
226
227 return 1;
228 }
229