1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
|
/*=========================================================================
Program: Visualization Toolkit
Module: vtkQuadraturePointsUtilities.hxx
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.
=========================================================================*/
#ifndef vtkQuadraturePointsUtilities_hxx
#define vtkQuadraturePointsUtilities_hxx
#include "vtkQuadratureSchemeDefinition.h"
#include "vtkDoubleArray.h"
#include "vtkFloatArray.h"
#include "vtkUnstructuredGrid.h"
namespace {
// Description:
// For all cells in the input "usg", for a specific array
// "V" interpolate to quadrature points using the given
// dictionary "dict" into "interpolated". Additionally if
// "indices" is not 0 then track the indices of where the
// values from each cell start as well. In the case of
// an error the return is 0.
template<class TV, class TI>
int Interpolate(
vtkUnstructuredGrid *usg,
const vtkIdType nCellsUsg,
TV *pV,
const int nCompsV,
vtkQuadratureSchemeDefinition **dict,
vtkDoubleArray *interpolated,
TI *indices)
{
// Walk cells.
vtkIdType currentIndex=0;
for (vtkIdType cellId=0; cellId<nCellsUsg; ++cellId)
{
// Point to the start of the data associated with this cell.
if (indices!=NULL)
{
indices[cellId]=static_cast<TI>(currentIndex);
}
// Grab the cell's associated shape function definition.
int cellType=usg->GetCellType(cellId);
vtkQuadratureSchemeDefinition *def=dict[cellType];
if (def==NULL)
{
// no quadrature scheme been specified for this cell type
// skipping the cell.
continue;
}
vtkIdType nNodes=def->GetNumberOfNodes();
int nQPts=def->GetNumberOfQuadraturePoints();
// Grab the cell's node ids.
vtkIdType *cellNodeIds=0;
usg->GetCellPoints(cellId,nNodes,cellNodeIds);
// Walk quadrature points.
for (int qPtId=0; qPtId<nQPts; ++qPtId)
{
// Grab the result and initialize.
double *r=interpolated->WritePointer(currentIndex,nCompsV);
for (int q=0; q<nCompsV; ++q)
{
r[q]=0.0;
}
// Grab shape function weights.
const double *N=def->GetShapeFunctionWeights(qPtId);
// Walk the cell's nodes.
for (int j=0; j<nNodes; ++j)
{
vtkIdType tupIdx=cellNodeIds[j]*nCompsV;
// Apply shape function weights.
for (int q=0; q<nCompsV; ++q)
{
r[q]+=N[j]*pV[tupIdx+q];
}
}
// Update the result index.
currentIndex+=nCompsV;
}
}
return 1;
}
// Description:
// Dispatch helper, descides what type of indices we are working with
template<class TV>
int Interpolate(
vtkUnstructuredGrid *usg,
const vtkIdType nCellsUsg,
TV *pV,
const int nCompsV,
vtkQuadratureSchemeDefinition **dict,
vtkDoubleArray *interpolated,
void *indices,
int indexType)
{
switch(indexType)
{
vtkTemplateMacro(
return Interpolate(usg,nCellsUsg,pV,nCompsV,dict,interpolated,(VTK_TT*)indices);
);
}
return 0;
}
//------------------------------------------------------------------------------
template <class T>
void ApplyShapeFunction(double *r,double N_j,T *A,int nComps)
{
for (int q=0; q<nComps; ++q)
{
r[q]+=N_j*A[q];
}
}
};
#endif
|
