1 /*=========================================================================
2 
3   Program:   Visualization Toolkit
4   Module:    vtkQuadratureSchemeDictionaryGenerator.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 
16 #include "vtkQuadratureSchemeDictionaryGenerator.h"
17 #include "vtkQuadratureSchemeDefinition.h"
18 
19 #include "vtkCellArray.h"
20 #include "vtkCellData.h"
21 #include "vtkCellTypes.h"
22 #include "vtkDataArray.h"
23 #include "vtkDoubleArray.h"
24 #include "vtkIdTypeArray.h"
25 #include "vtkInformation.h"
26 #include "vtkInformationQuadratureSchemeDefinitionVectorKey.h"
27 #include "vtkInformationVector.h"
28 #include "vtkIntArray.h"
29 #include "vtkObjectFactory.h"
30 #include "vtkPointData.h"
31 #include "vtkPoints.h"
32 #include "vtkPolyData.h"
33 #include "vtkType.h"
34 #include "vtkUnstructuredGrid.h"
35 #include "vtkUnstructuredGridAlgorithm.h"
36 
37 #include "vtkSmartPointer.h"
38 #include <sstream>
39 #include <string>
40 using std::ostringstream;
41 using std::string;
42 
43 // Here are some default shape functions weights which
44 // we will use to create dictionaries in a gvien data set.
45 // Unused weights are commented out to avoid compiler warnings.
46 namespace
47 {
48 // double W_T_11_A[]={
49 //     3.33333333333334e-01, 3.33333333333333e-01, 3.33333333333333e-01};
50 
51 double W_T_32_A[] = { 1.66666666666660e-01, 6.66666666666670e-01, 1.66666666666670e-01,
52   6.66666666666660e-01, 1.66666666666670e-01, 1.66666666666670e-01, 1.66666666666660e-01,
53   1.66666666666670e-01, 6.66666666666670e-01 };
54 
55 // double W_T_32_B[]={
56 //     5.00000000000000e-01, 5.00000000000000e-01, 0.00000000000000e+00,
57 //     5.00000000000000e-01, 0.00000000000000e+00, 5.00000000000000e-01,
58 //     0.00000000000000e+00, 5.00000000000000e-01, 5.00000000000000e-01};
59 
60 double W_QT_43_A[] = { -1.11111111111111e-01, -1.11111111111111e-01, -1.11111111111111e-01,
61   4.44444444444445e-01, 4.44444444444444e-01, 4.44444444444445e-01, -1.20000000000000e-01,
62   1.20000000000000e-01, -1.20000000000000e-01, 4.80000000000000e-01, 4.80000000000000e-01,
63   1.60000000000000e-01, 1.20000000000000e-01, -1.20000000000000e-01, -1.20000000000000e-01,
64   4.80000000000000e-01, 1.60000000000000e-01, 4.80000000000000e-01, -1.20000000000000e-01,
65   -1.20000000000000e-01, 1.20000000000000e-01, 1.60000000000000e-01, 4.80000000000000e-01,
66   4.80000000000000e-01 };
67 
68 double W_Q_42_A[] = { 6.22008467928145e-01, 1.66666666666667e-01, 4.46581987385206e-02,
69   1.66666666666667e-01, 1.66666666666667e-01, 4.46581987385206e-02, 1.66666666666667e-01,
70   6.22008467928145e-01, 1.66666666666667e-01, 6.22008467928145e-01, 1.66666666666667e-01,
71   4.46581987385206e-02, 4.46581987385206e-02, 1.66666666666667e-01, 6.22008467928145e-01,
72   1.66666666666667e-01 };
73 
74 double W_QQ_93_A[] = { 4.32379000772438e-01, -1.00000000000001e-01, -3.23790007724459e-02,
75   -1.00000000000001e-01, 3.54919333848301e-01, 4.50806661517046e-02, 4.50806661517046e-02,
76   3.54919333848301e-01, -1.00000000000001e-01, -1.00000000000001e-01, -1.00000000000001e-01,
77   -1.00000000000001e-01, 2.00000000000003e-01, 1.12701665379260e-01, 2.00000000000003e-01,
78   8.87298334620740e-01, -1.00000000000001e-01, -3.23790007724459e-02, -1.00000000000001e-01,
79   4.32379000772438e-01, 4.50806661517046e-02, 4.50806661517046e-02, 3.54919333848301e-01,
80   3.54919333848301e-01, -1.00000000000001e-01, -1.00000000000001e-01, -1.00000000000001e-01,
81   -1.00000000000001e-01, 8.87298334620740e-01, 2.00000000000003e-01, 1.12701665379260e-01,
82   2.00000000000003e-01, -2.50000000000000e-01, -2.50000000000000e-01, -2.50000000000000e-01,
83   -2.50000000000000e-01, 5.00000000000000e-01, 5.00000000000000e-01, 5.00000000000000e-01,
84   5.00000000000000e-01, -1.00000000000001e-01, -1.00000000000001e-01, -1.00000000000001e-01,
85   -1.00000000000001e-01, 1.12701665379260e-01, 2.00000000000003e-01, 8.87298334620740e-01,
86   2.00000000000003e-01, -1.00000000000001e-01, 4.32379000772438e-01, -1.00000000000001e-01,
87   -3.23790007724459e-02, 3.54919333848301e-01, 3.54919333848301e-01, 4.50806661517046e-02,
88   4.50806661517046e-02, -1.00000000000001e-01, -1.00000000000001e-01, -1.00000000000001e-01,
89   -1.00000000000001e-01, 2.00000000000003e-01, 8.87298334620740e-01, 2.00000000000003e-01,
90   1.12701665379260e-01, -3.23790007724459e-02, -1.00000000000001e-01, 4.32379000772438e-01,
91   -1.00000000000001e-01, 4.50806661517046e-02, 3.54919333848301e-01, 3.54919333848301e-01,
92   4.50806661517046e-02 };
93 
94 // double W_E_41_A[]={
95 //      2.50000000000000e-01, 2.50000000000000e-01, 2.50000000000000e-01, 2.50000000000000e-01};
96 
97 double W_E_42_A[] = { 6.25000000000000e-01, 1.25000000000000e-01, 1.25000000000000e-01,
98   1.25000000000000e-01, 1.25000000000000e-01, 5.62500000000000e-01, 1.87500000000000e-01,
99   1.25000000000000e-01, 1.25000000000000e-01, 1.87500000000000e-01, 5.62500000000000e-01,
100   1.25000000000000e-01, 1.25000000000000e-01, 6.25000000000000e-02, 6.25000000000000e-02,
101   7.50000000000000e-01 };
102 
103 // double W_QE_41_A[]={
104 //     -1.25000000000000e-01, -1.25000000000000e-01, -1.25000000000000e-01,
105 //     -1.25000000000000e-01, 2.50000000000000e-01, 2.50000000000000e-01, 2.50000000000000e-01, 2.50000000000000e-01,
106 //     2.50000000000000e-01, 2.50000000000000e-01};
107 
108 double W_QE_42_A[] = { 1.56250000000000e-01, -9.37500000000000e-02, -9.37500000000000e-02,
109   -9.37500000000000e-02, 3.12500000000000e-01, 6.25000000000000e-02, 3.12500000000000e-01,
110   3.12500000000000e-01, 6.25000000000000e-02, 6.25000000000000e-02, -9.37500000000000e-02,
111   7.03125000000000e-02, -1.17187500000000e-01, -9.37500000000000e-02, 2.81250000000000e-01,
112   4.21875000000000e-01, 9.37500000000000e-02, 6.25000000000000e-02, 2.81250000000000e-01,
113   9.37500000000000e-02, -9.37500000000000e-02, -1.17187500000000e-01, 7.03125000000000e-02,
114   -9.37500000000000e-02, 9.37500000000000e-02, 4.21875000000000e-01, 2.81250000000000e-01,
115   6.25000000000000e-02, 9.37500000000000e-02, 2.81250000000000e-01, -9.37500000000000e-02,
116   -5.46875000000000e-02, -5.46875000000000e-02, 3.75000000000000e-01, 3.12500000000000e-02,
117   1.56250000000000e-02, 3.12500000000000e-02, 3.75000000000000e-01, 1.87500000000000e-01,
118   1.87500000000000e-01 };
119 
120 };
121 
122 vtkStandardNewMacro(vtkQuadratureSchemeDictionaryGenerator);
123 
124 //------------------------------------------------------------------------------
vtkQuadratureSchemeDictionaryGenerator()125 vtkQuadratureSchemeDictionaryGenerator::vtkQuadratureSchemeDictionaryGenerator()
126 {
127   this->SetNumberOfInputPorts(1);
128   this->SetNumberOfOutputPorts(1);
129 }
130 
131 //------------------------------------------------------------------------------
132 vtkQuadratureSchemeDictionaryGenerator::~vtkQuadratureSchemeDictionaryGenerator() = default;
133 
134 //------------------------------------------------------------------------------
FillInputPortInformation(int port,vtkInformation * info)135 int vtkQuadratureSchemeDictionaryGenerator::FillInputPortInformation(int port, vtkInformation* info)
136 {
137   switch (port)
138   {
139     case 0:
140       info->Set(vtkDataObject::DATA_TYPE_NAME(), "vtkUnstructuredGrid");
141       break;
142   }
143   return 1;
144 }
145 
146 //------------------------------------------------------------------------------
FillOutputPortInformation(int port,vtkInformation * info)147 int vtkQuadratureSchemeDictionaryGenerator::FillOutputPortInformation(
148   int port, vtkInformation* info)
149 {
150   switch (port)
151   {
152     case 0:
153       info->Set(vtkDataObject::DATA_TYPE_NAME(), "vtkUnstructuredGrid");
154       break;
155   }
156   return 1;
157 }
158 
159 //------------------------------------------------------------------------------
RequestData(vtkInformation *,vtkInformationVector ** input,vtkInformationVector * output)160 int vtkQuadratureSchemeDictionaryGenerator::RequestData(
161   vtkInformation*, vtkInformationVector** input, vtkInformationVector* output)
162 {
163   vtkDataObject* tmpDataObj;
164   // Get the inputs
165   tmpDataObj = input[0]->GetInformationObject(0)->Get(vtkDataObject::DATA_OBJECT());
166   vtkUnstructuredGrid* usgIn = vtkUnstructuredGrid::SafeDownCast(tmpDataObj);
167   // Get the outputs
168   tmpDataObj = output->GetInformationObject(0)->Get(vtkDataObject::DATA_OBJECT());
169   vtkUnstructuredGrid* usgOut = vtkUnstructuredGrid::SafeDownCast(tmpDataObj);
170 
171   // Quick sanity check.
172   if (usgIn == nullptr || usgOut == nullptr || usgIn->GetNumberOfPoints() == 0 ||
173     usgIn->GetPointData()->GetNumberOfArrays() == 0)
174   {
175     vtkWarningMacro("Filter data has not been configured correctly. Aborting.");
176     return 1;
177   }
178 
179   // Copy the unstructured grid on the input
180   usgOut->ShallowCopy(usgIn);
181 
182   // Interpolate the data arrays, but no points. Results
183   // are stored in field data arrays.
184   this->Generate(usgOut);
185 
186   return 1;
187 }
188 
189 //------------------------------------------------------------------------------
Generate(vtkUnstructuredGrid * usgOut)190 int vtkQuadratureSchemeDictionaryGenerator::Generate(vtkUnstructuredGrid* usgOut)
191 {
192   // Get the dictionary key.
193   vtkInformationQuadratureSchemeDefinitionVectorKey* key =
194     vtkQuadratureSchemeDefinition::DICTIONARY();
195 
196   // Get the cell types used by the data set.
197   vtkCellTypes* cellTypes = vtkCellTypes::New();
198   usgOut->GetCellTypes(cellTypes);
199   // add a definition to the dictionary for each cell type.
200   int nCellTypes = cellTypes->GetNumberOfTypes();
201 
202   // create the offset array and store the dictionary within
203   vtkIdTypeArray* offsets = vtkIdTypeArray::New();
204   string basename = "QuadratureOffset";
205   string finalname = basename;
206   vtkDataArray* data = usgOut->GetCellData()->GetArray(basename.c_str());
207   ostringstream interpolatedName;
208   int i = 0;
209   while (data != nullptr)
210   {
211     interpolatedName << basename << i;
212     data = usgOut->GetCellData()->GetArray(interpolatedName.str().c_str());
213     finalname = interpolatedName.str();
214     i++;
215   }
216 
217   offsets->SetName(finalname.c_str());
218   usgOut->GetCellData()->AddArray(offsets);
219   vtkInformation* info = offsets->GetInformation();
220 
221   for (int typeId = 0; typeId < nCellTypes; ++typeId)
222   {
223     int cellType = cellTypes->GetCellType(typeId);
224     // Initiaze a definition for this particular cell type.
225     vtkSmartPointer<vtkQuadratureSchemeDefinition> def =
226       vtkSmartPointer<vtkQuadratureSchemeDefinition>::New();
227     switch (cellType)
228     {
229       case VTK_TRIANGLE:
230         def->Initialize(VTK_TRIANGLE, 3, 3, W_T_32_A);
231         break;
232       case VTK_QUADRATIC_TRIANGLE:
233         def->Initialize(VTK_QUADRATIC_TRIANGLE, 6, 4, W_QT_43_A);
234         break;
235       case VTK_QUAD:
236         def->Initialize(VTK_QUAD, 4, 4, W_Q_42_A);
237         break;
238       case VTK_QUADRATIC_QUAD:
239         def->Initialize(VTK_QUADRATIC_QUAD, 8, 9, W_QQ_93_A);
240         break;
241       case VTK_TETRA:
242         def->Initialize(VTK_TETRA, 4, 4, W_E_42_A);
243         break;
244       case VTK_QUADRATIC_TETRA:
245         def->Initialize(VTK_QUADRATIC_TETRA, 10, 4, W_QE_42_A);
246         break;
247       default:
248         cerr << "Error: Cell type " << cellType << " found "
249              << "with no definition provided. Add a definition "
250              << " in " << __FILE__ << ". Aborting." << endl;
251         return 0;
252     }
253 
254     // The definition must appear in the dictionary associated
255     // with the offset array
256     key->Set(info, def, cellType);
257   }
258 
259   int dictSize = key->Size(info);
260   vtkQuadratureSchemeDefinition** dict = new vtkQuadratureSchemeDefinition*[dictSize];
261   key->GetRange(info, dict, 0, 0, dictSize);
262 
263   offsets->SetNumberOfTuples(usgOut->GetNumberOfCells());
264   vtkIdType offset = 0;
265   for (vtkIdType cellid = 0; cellid < usgOut->GetNumberOfCells(); cellid++)
266   {
267     offsets->SetValue(cellid, offset);
268     vtkCell* cell = usgOut->GetCell(cellid);
269     int cellType = cell->GetCellType();
270     vtkQuadratureSchemeDefinition* celldef = dict[cellType];
271     offset += celldef->GetNumberOfQuadraturePoints();
272   }
273   offsets->Delete();
274   cellTypes->Delete();
275   delete[] dict;
276   return 1;
277 }
278 
279 //------------------------------------------------------------------------------
PrintSelf(ostream & os,vtkIndent indent)280 void vtkQuadratureSchemeDictionaryGenerator::PrintSelf(ostream& os, vtkIndent indent)
281 {
282   this->Superclass::PrintSelf(os, indent);
283 
284   os << indent << "No state." << endl;
285 }
286