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