1 #include "orthogonalityoptimiser.h"
2 #include "surfaceprojection.h"
3 #include "guimainwindow.h"
5 OrthogonalityOptimiser::OrthogonalityOptimiser()
7 m_StrictPrismChecking
= true;
10 double OrthogonalityOptimiser::faceError(const face_t
&face
, vec3_t x
)
12 vec3_t n
= face
.c1
- face
.c2
;
15 foreach (vtkIdType id_node
, face
.nodes
) {
16 if (id_node
!= m_IdNode
) {
18 m_Grid
->GetPoint(id_node
, xn
.data());
22 xm
*= 1.0/face
.nodes
.size();
24 if (!checkVector(v
)) {
30 double OrthogonalityOptimiser::func(vec3_t x
)
32 if (m_IdNode
== 427) {
33 //cout << "stop" << endl;
36 foreach (int i_face
, m_Node2Face
[m_IdNode
]) {
37 double e
= faceError(m_Faces
[i_face
], x
);
49 void OrthogonalityOptimiser::getFace(vtkIdType id_cell
, int subidx
, QSet
<vtkIdType
> &nodes
)
52 vtkIdType N_pts
, *pts
;
53 m_Grid
->GetCellPoints(id_cell
, N_pts
, pts
);
54 vtkIdType type_cell
= m_Grid
->GetCellType(id_cell
);
55 if (type_cell
== VTK_TETRA
) {
77 if (type_cell
== VTK_WEDGE
) {
101 nodes
.insert(pts
[0]);
102 nodes
.insert(pts
[2]);
103 nodes
.insert(pts
[5]);
104 nodes
.insert(pts
[3]);
107 if (type_cell
== VTK_HEXAHEDRON
) {
109 nodes
.insert(pts
[0]);
110 nodes
.insert(pts
[3]);
111 nodes
.insert(pts
[2]);
112 nodes
.insert(pts
[1]);
115 nodes
.insert(pts
[4]);
116 nodes
.insert(pts
[5]);
117 nodes
.insert(pts
[6]);
118 nodes
.insert(pts
[7]);
121 nodes
.insert(pts
[0]);
122 nodes
.insert(pts
[1]);
123 nodes
.insert(pts
[5]);
124 nodes
.insert(pts
[4]);
127 nodes
.insert(pts
[3]);
128 nodes
.insert(pts
[7]);
129 nodes
.insert(pts
[6]);
130 nodes
.insert(pts
[2]);
133 nodes
.insert(pts
[0]);
134 nodes
.insert(pts
[4]);
135 nodes
.insert(pts
[7]);
136 nodes
.insert(pts
[3]);
139 nodes
.insert(pts
[1]);
140 nodes
.insert(pts
[2]);
141 nodes
.insert(pts
[6]);
142 nodes
.insert(pts
[5]);
147 int OrthogonalityOptimiser::getNumberOfFaces(vtkIdType id_cell
)
149 vtkIdType cell_type
= m_Grid
->GetCellType(id_cell
);
151 if (cell_type
== VTK_TETRA
) N
= 4;
152 else if (cell_type
== VTK_PYRAMID
) N
= 5;
153 else if (cell_type
== VTK_WEDGE
) N
= 5;
154 else if (cell_type
== VTK_HEXAHEDRON
) N
= 6;
158 void OrthogonalityOptimiser::buildNode2Face()
160 l2g_t cells
= m_Part
.getCells();
161 m_Node2Face
.resize(m_Grid
->GetNumberOfPoints());
163 foreach (vtkIdType id_cell1
, cells
) {
164 if (!isSurface(id_cell1
, m_Grid
)) {
165 for (int i
= 0; i
< m_Part
.c2cGSize(id_cell1
); ++i
) {
166 vtkIdType id_cell2
= m_Part
.c2cGG(id_cell1
, i
);
167 if (id_cell1
< id_cell2
&& !isSurface(id_cell2
, m_Grid
)) {
168 for (int j
= 0; j
< getNumberOfFaces(id_cell1
); ++j
) {
169 QSet
<vtkIdType
> nodes1
;
170 getFace(id_cell1
, j
, nodes1
);
171 for (int k
= 0; k
< getNumberOfFaces(id_cell2
); ++k
) {
172 QSet
<vtkIdType
> nodes2
;
173 getFace(id_cell2
, j
, nodes2
);
174 if (nodes1
== nodes2
) {
176 F
.c1
= cellCentre(m_Grid
, id_cell1
);
177 F
.c2
= cellCentre(m_Grid
, id_cell2
);
178 F
.nodes
.resize(nodes1
.size());
179 qCopy(nodes1
.begin(), nodes1
.end(), F
.nodes
.begin());
188 m_Faces
.resize(faces
.size());
189 qCopy(faces
.begin(), faces
.end(), m_Faces
.begin());
190 for (int i
= 0; i
< m_Faces
.size(); ++i
) {
191 foreach (vtkIdType id_node
, m_Faces
[i
].nodes
) {
192 m_Node2Face
[id_node
].append(i
);
197 void OrthogonalityOptimiser::fixBoundaryNodes()
199 l2g_t cells
= m_Part
.getCells();
200 m_Fixed
.fill(false, m_Grid
->GetNumberOfPoints());
201 EG_VTKDCC(vtkIntArray
, cell_code
, m_Grid
, "cell_code");
202 QVector
<QSet
<int> > n2bc(m_Grid
->GetNumberOfPoints());
203 foreach (vtkIdType id_cell
, cells
) {
204 if (isSurface(id_cell
, m_Grid
)) {
205 vtkIdType N_pts
, *pts
;
206 m_Grid
->GetCellPoints(id_cell
, N_pts
, pts
);
207 for (int i
= 0; i
< N_pts
; ++i
) {
208 n2bc
[pts
[i
]].insert(cell_code
->GetValue(id_cell
));
212 m_Projection
.fill(NULL
, m_Grid
->GetNumberOfPoints());
213 for (vtkIdType id_node
= 0; id_node
< m_Grid
->GetNumberOfPoints(); ++id_node
) {
214 if (n2bc
[id_node
].size() > 1) {
215 m_Fixed
[id_node
] = true;
216 } else if (n2bc
[id_node
].size() == 1) {
217 GuiMainWindow::pointer()->getSurfProj(n2bc
[id_node
].toList().first());
219 if (m_Node2Face
[id_node
].size() == 0) {
220 m_Fixed
[id_node
] = true;
225 vec3_t
OrthogonalityOptimiser::newPosition(vtkIdType id_node
)
228 m_Grid
->GetPoint(id_node
, x_old
.data());
229 if (m_Node2Face
[id_node
].size() == 0) {
233 foreach (int i_face
, m_Node2Face
[id_node
]) {
234 const face_t
&face
= m_Faces
[i_face
];
235 vec3_t n
= face
.c1
- face
.c2
;
238 foreach (vtkIdType id_node
, face
.nodes
) {
240 m_Grid
->GetPoint(id_node
, xn
.data());
243 xm
*= 1.0/face
.nodes
.size();
244 vec3_t v
= x_old
- xm
;
245 if (!checkVector(v
)) {
251 x_new
*= (1.0/m_Node2Face
[id_node
].size());
255 void OrthogonalityOptimiser::operate()
260 l2g_t nodes
= m_Part
.getNodes();
262 foreach (vtkIdType id_node
, nodes
) {
263 if (!m_Fixed
[id_node
]) {
266 m_Grid
->GetPoint(id_node
, x_old
.data());
267 //vec3_t x_new = optimise(x_old);
268 vec3_t x_new
= newPosition(id_node
);
269 if (checkVector(x_new
)) {
270 if (m_Projection
[id_node
]) {
271 x_new
= m_Projection
[id_node
]->projectRestricted(x_new
, id_node
);
273 if (checkVector(x_new
)) {
274 vec3_t Dx
= x_new
- x_old
;
275 if (moveNode(id_node
, Dx
)) ++N
;
280 cout
<< N
<< " nodes optimised" << endl
;