2 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
4 // + This file is part of enGrid. +
6 // + Copyright 2008,2009 Oliver Gloth +
8 // + enGrid is free software: you can redistribute it and/or modify +
9 // + it under the terms of the GNU General Public License as published by +
10 // + the Free Software Foundation, either version 3 of the License, or +
11 // + (at your option) any later version. +
13 // + enGrid is distributed in the hope that it will be useful, +
14 // + but WITHOUT ANY WARRANTY; without even the implied warranty of +
15 // + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +
16 // + GNU General Public License for more details. +
18 // + You should have received a copy of the GNU General Public License +
19 // + along with enGrid. If not, see <http://www.gnu.org/licenses/>. +
21 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
23 #include "operation.h"
24 #include "guimainwindow.h"
26 #include "geometrytools.h"
27 using namespace GeometryTools
;
29 #include <QApplication>
31 QSet
<Operation
*> Operation::garbage_operations
;
33 void Operation::collectGarbage()
35 QSet
<Operation
*> delete_operations
;
36 foreach (Operation
*op
, garbage_operations
) {
37 if (!op
->getThread().isRunning()) {
38 delete_operations
.insert(op
);
39 cout
<< "deleting Operation " << op
<< endl
;
43 foreach (Operation
*op
, delete_operations
) {
44 garbage_operations
.remove(op
);
48 Operation::Operation()
56 Operation::~Operation()
66 garbage_operations
.insert(this);
69 void OperationThread::run()
72 GuiMainWindow::lock();
73 GuiMainWindow::pointer()->setBusy();
76 op
->err
= new Error();
79 GuiMainWindow::unlock();
80 GuiMainWindow::pointer()->setIdle();
83 void Operation::operator()()
86 if (GuiMainWindow::tryLock()) {
88 thread
.setOperation(this);
89 GuiMainWindow::unlock();
90 thread
.start(QThread::LowPriority
);
92 QMessageBox::warning(NULL
, "not permitted", "Operation is not permitted while background process is running!");
104 void Operation::setAllCells()
106 QVector
<vtkIdType
> all_cells
;
107 getAllCells(all_cells
, grid
);
111 void Operation::setAllVolumeCells()
113 QVector
<vtkIdType
> cells
;
114 getAllVolumeCells(cells
, grid
);
118 void Operation::setAllSurfaceCells()
120 QVector
<vtkIdType
> cells
;
121 getAllSurfaceCells(cells
, grid
);
125 void Operation::initMapping()
127 nodes_map
.resize(nodes
.size());
128 for (int i_nodes
= 0; i_nodes
< nodes
.size(); ++i_nodes
) {
129 nodes_map
[i_nodes
] = nodes
[i_nodes
];
131 cells_map
.resize(cells
.size());
132 for (int i_cells
= 0; i_cells
< cells
.size(); ++i_cells
) {
133 cells_map
[i_cells
] = cells
[i_cells
];
137 void Operation::checkGrid()
140 grid
= GuiMainWindow::pointer()->getGrid();
142 if ((cells
.size() == 0) && autoset
) {
147 void Operation::updateActors()
149 mainWindow()->updateActors();
152 GuiMainWindow
* Operation::mainWindow()
154 return GuiMainWindow::pointer();
157 void Operation::populateBoundaryCodes(QListWidget
*lw
)
160 mainWindow()->getAllBoundaryCodes(bcs
);
161 foreach(int bc
, bcs
) {
162 QListWidgetItem
*lwi
= new QListWidgetItem(lw
);
163 lwi
->setCheckState(Qt::Unchecked
);
165 QTextStream
ts(&text
);
168 lwi
->setFlags(Qt::ItemIsUserCheckable
| Qt::ItemIsEnabled
);
172 stencil_t
Operation::getStencil(vtkIdType id_cell1
, int j1
)
176 S
.id_cell1
= id_cell1
;
177 if (c2c
[_cells
[id_cell1
]][j1
] != -1) {
178 S
.id_cell2
= cells
[c2c
[_cells
[id_cell1
]][j1
]];
179 if (grid
->GetCellType(S
.id_cell2
) != VTK_TRIANGLE
) {
182 vtkIdType N1
, N2
, *pts1
, *pts2
;
183 grid
->GetCellPoints(S
.id_cell1
, N1
, pts1
);
184 grid
->GetCellPoints(S
.id_cell2
, N2
, pts2
);
185 if (j1
== 0) { S
.p
[0] = pts1
[2]; S
.p
[1] = pts1
[0]; S
.p
[3] = pts1
[1]; }
186 else if (j1
== 1) { S
.p
[0] = pts1
[0]; S
.p
[1] = pts1
[1]; S
.p
[3] = pts1
[2]; }
187 else if (j1
== 2) { S
.p
[0] = pts1
[1]; S
.p
[1] = pts1
[2]; S
.p
[3] = pts1
[0]; };
189 if (c2c
[_cells
[S
.id_cell2
]][0] != -1) {
190 if (cells
[c2c
[_cells
[S
.id_cell2
]][0]] == S
.id_cell1
) {
195 if (c2c
[_cells
[S
.id_cell2
]][1] != -1) {
196 if (cells
[c2c
[_cells
[S
.id_cell2
]][1]] == S
.id_cell1
) {
201 if (c2c
[_cells
[S
.id_cell2
]][2] != -1) {
202 if (cells
[c2c
[_cells
[S
.id_cell2
]][2]] == S
.id_cell1
) {
214 grid
->GetCellPoints(S
.id_cell1
, N1
, pts1
);
215 if (j1
== 0) { S
.p
[0] = pts1
[2]; S
.p
[1] = pts1
[0]; S
.p
[3] = pts1
[1]; }
216 else if (j1
== 1) { S
.p
[0] = pts1
[0]; S
.p
[1] = pts1
[1]; S
.p
[3] = pts1
[2]; }
217 else if (j1
== 2) { S
.p
[0] = pts1
[1]; S
.p
[1] = pts1
[2]; S
.p
[3] = pts1
[0]; };
222 ostream
& operator<<(ostream
&out
, stencil_t S
)
224 out
<<"S.id_cell1="<<S
.id_cell1
<<" ";
225 out
<<"S.id_cell2="<<S
.id_cell2
<<" ";
226 out
<<"S.valid="<<S
.valid
<<" ";
228 for(int i
=0;i
<4;i
++){
236 //////////////////////////////////////////////
237 double CurrentVertexAvgDist(vtkIdType a_vertex
,QVector
< QSet
< int > > &n2n
,vtkUnstructuredGrid
*a_grid
)
241 int N
=n2n
[a_vertex
].size();
243 a_grid
->GetPoint(a_vertex
, C
.data());
244 foreach(int i
,n2n
[a_vertex
])
247 a_grid
->GetPoint(i
, M
.data());
248 total_dist
+=(M
-C
).abs();
250 avg_dist
=total_dist
/(double)N
;
254 double CurrentMeshDensity(vtkIdType a_vertex
,QVector
< QSet
< int > > &n2n
,vtkUnstructuredGrid
*a_grid
)
258 int N
=n2n
[a_vertex
].size();
260 a_grid
->GetPoint(a_vertex
, C
.data());
261 foreach(int i
,n2n
[a_vertex
])
264 a_grid
->GetPoint(i
, M
.data());
265 total_dist
+=(M
-C
).abs();
267 avg_dist
=total_dist
/(double)N
;
268 double avg_density
=1./avg_dist
;
272 double DesiredVertexAvgDist(vtkIdType a_vertex
,QVector
< QSet
< int > > &n2n
,vtkUnstructuredGrid
*a_grid
)
276 EG_VTKDCN(vtkDoubleArray
, node_meshdensity
, a_grid
, "node_meshdensity");
277 int N
=n2n
[a_vertex
].size();
278 foreach(int i
,n2n
[a_vertex
])
280 total_dist
+=1./node_meshdensity
->GetValue(i
);
282 avg_dist
=total_dist
/(double)N
;
286 double DesiredMeshDensity(vtkIdType a_vertex
,QVector
< QSet
< int > > &n2n
,vtkUnstructuredGrid
*a_grid
)
288 double total_density
=0;
289 double avg_density
=0;
290 EG_VTKDCN(vtkDoubleArray
, node_meshdensity
, a_grid
, "node_meshdensity");
291 int N
=n2n
[a_vertex
].size();
292 foreach(int i
,n2n
[a_vertex
])
294 total_density
+=node_meshdensity
->GetValue(i
);
296 avg_density
=total_density
/(double)N
;
300 ///////////////////////////////////////////
301 vtkIdType
Operation::getClosestNode(vtkIdType a_id_node
,vtkUnstructuredGrid
* a_grid
)
304 a_grid
->GetPoint(a_id_node
,C
.data());
305 vtkIdType id_minlen
=-1;
307 foreach(vtkIdType neighbour
,n2n
[a_id_node
])
310 a_grid
->GetPoint(neighbour
,M
.data());
311 double len
=(M
-C
).abs();
312 if(minlen
<0 or len
<minlen
)
321 vtkIdType
Operation::getFarthestNode(vtkIdType a_id_node
,vtkUnstructuredGrid
* a_grid
)
324 a_grid
->GetPoint(a_id_node
,C
.data());
325 vtkIdType id_maxlen
=-1;
327 foreach(vtkIdType neighbour
,n2n
[a_id_node
])
330 a_grid
->GetPoint(neighbour
,M
.data());
331 double len
=(M
-C
).abs();
332 if(maxlen
<0 or len
>maxlen
)
341 bool Operation::SwapCells(vtkUnstructuredGrid
* a_grid
, stencil_t S
)
345 vec3_t x3
[4], x3_0(0,0,0);
347 for (int k
= 0; k
< 4; ++k
) {
348 a_grid
->GetPoints()->GetPoint(S
.p
[k
], x3
[k
].data());
351 vec3_t n1
= triNormal(x3
[0], x3
[1], x3
[3]);
352 vec3_t n2
= triNormal(x3
[1], x3
[2], x3
[3]);
355 if ( (n1
*n2
) > 0.8) {
358 vec3_t ex
= orthogonalVector(n
);
359 vec3_t ey
= ex
.cross(n
);
360 for (int k
= 0; k
< 4; ++k
) {
361 x
[k
] = vec2_t(x3
[k
]*ex
, x3
[k
]*ey
);
363 vec2_t r1
, r2
, r3
, u1
, u2
, u3
;
364 r1
= 0.5*(x
[0] + x
[1]); u1
= turnLeft(x
[1] - x
[0]);
365 r2
= 0.5*(x
[1] + x
[2]); u2
= turnLeft(x
[2] - x
[1]);
366 r3
= 0.5*(x
[1] + x
[3]); u3
= turnLeft(x
[3] - x
[1]);
370 if (intersection(k
, l
, r1
, u1
, r3
, u3
)) {
372 if (intersection(k
, l
, r2
, u2
, r3
, u3
)) {
382 if ((xm1
- x
[2]).abs() < (xm1
- x
[0]).abs()) {
385 if ((xm2
- x
[0]).abs() < (xm2
- x
[2]).abs()) {
392 vtkIdType new_pts1
[3], new_pts2
[3];
393 new_pts1
[0] = S
.p
[1];
394 new_pts1
[1] = S
.p
[2];
395 new_pts1
[2] = S
.p
[0];
396 new_pts2
[0] = S
.p
[2];
397 new_pts2
[1] = S
.p
[3];
398 new_pts2
[2] = S
.p
[0];
399 a_grid
->ReplaceCell(S
.id_cell1
, 3, new_pts1
);
400 a_grid
->ReplaceCell(S
.id_cell2
, 3, new_pts2
);
405 void Operation::quad2triangle(vtkUnstructuredGrid
* src
,vtkIdType quadcell
)
407 vtkIdType type_cell
= grid
->GetCellType(quadcell
);
408 cout
<<"It's a "<<type_cell
<<endl
;
409 if(type_cell
==VTK_QUAD
)
411 cout_grid(cout
,src
,true,true,true,true);
412 EG_VTKSP(vtkUnstructuredGrid
, dst
);
414 int N_points
=src
->GetNumberOfPoints();
415 int N_cells
=src
->GetNumberOfCells();
416 allocateGrid(dst
,N_cells
+1,N_points
);
418 for (vtkIdType id_node
= 0; id_node
< src
->GetNumberOfPoints(); ++id_node
) {
420 src
->GetPoints()->GetPoint(id_node
, x
.data());
421 dst
->GetPoints()->SetPoint(id_node
, x
.data());
422 copyNodeData(src
, id_node
, dst
, id_node
);
424 for (vtkIdType id_cell
= 0; id_cell
< src
->GetNumberOfCells(); ++id_cell
) {
425 vtkIdType N_pts
, *pts
;
426 vtkIdType type_cell
= src
->GetCellType(id_cell
);
427 src
->GetCellPoints(id_cell
, N_pts
, pts
);
428 vtkIdType id_new_cell
;
429 vtkIdType id_new_cell1
;
430 vtkIdType id_new_cell2
;
431 if(id_cell
!=quadcell
)
433 id_new_cell
= dst
->InsertNextCell(type_cell
, N_pts
, pts
);
434 copyCellData(src
, id_cell
, dst
, id_new_cell
);
438 vtkIdType triangle1
[3];
439 vtkIdType triangle2
[3];
446 id_new_cell1
= dst
->InsertNextCell(VTK_TRIANGLE
, 3, triangle1
);
447 copyCellData(src
, id_cell
, dst
, id_new_cell1
);
448 id_new_cell2
= dst
->InsertNextCell(VTK_TRIANGLE
, 3, triangle2
);
449 copyCellData(src
, id_cell
, dst
, id_new_cell2
);
451 S
.id_cell1
=id_new_cell1
;
452 S
.id_cell2
=id_new_cell2
;
461 cout_grid(cout
,dst
,true,true,true,true);
466 void Operation::quad2triangle(vtkUnstructuredGrid
* src
,vtkIdType quadcell
,vtkIdType MovingPoint
)
468 vtkIdType type_cell
= grid
->GetCellType(quadcell
);
469 cout
<<"It's a "<<type_cell
<<endl
;
470 if(type_cell
==VTK_QUAD
)
472 cout_grid(cout
,src
,true,true,true,true);
473 EG_VTKSP(vtkUnstructuredGrid
, dst
);
475 int N_points
=src
->GetNumberOfPoints();
476 int N_cells
=src
->GetNumberOfCells();
477 allocateGrid(dst
,N_cells
+1,N_points
);
479 for (vtkIdType id_node
= 0; id_node
< src
->GetNumberOfPoints(); ++id_node
) {
481 src
->GetPoints()->GetPoint(id_node
, x
.data());
482 dst
->GetPoints()->SetPoint(id_node
, x
.data());
483 copyNodeData(src
, id_node
, dst
, id_node
);
485 for (vtkIdType id_cell
= 0; id_cell
< src
->GetNumberOfCells(); ++id_cell
) {
486 vtkIdType N_pts
, *pts
;
487 src
->GetCellPoints(id_cell
, N_pts
, pts
);
488 vtkIdType type_cell
= src
->GetCellType(id_cell
);
489 vtkIdType id_new_cell
;
490 vtkIdType id_new_cell1
;
491 vtkIdType id_new_cell2
;
492 if(id_cell
!=quadcell
)
494 id_new_cell
= dst
->InsertNextCell(type_cell
, N_pts
, pts
);
495 copyCellData(src
, id_cell
, dst
, id_new_cell
);
499 vtkIdType triangle1
[3];
500 vtkIdType triangle2
[3];
501 if(MovingPoint
==pts
[0] || MovingPoint
==pts
[2])
519 id_new_cell1
= dst
->InsertNextCell(VTK_TRIANGLE
, 3, triangle1
);
520 copyCellData(src
, id_cell
, dst
, id_new_cell1
);
521 id_new_cell2
= dst
->InsertNextCell(VTK_TRIANGLE
, 3, triangle2
);
522 copyCellData(src
, id_cell
, dst
, id_new_cell2
);
525 cout_grid(cout
,dst
,true,true,true,true);
530 int Operation::NumberOfCommonPoints(vtkIdType node1
, vtkIdType node2
, bool& IsTetra
)
532 // QVector< QSet< int > > n2n
533 QSet
<int> node1_neighbours
=n2n
[node1
];
534 QSet
<int> node2_neighbours
=n2n
[node2
];
535 QSet
<int> intersection
=node1_neighbours
.intersect(node2_neighbours
);
536 int N
=intersection
.size();
540 QSet
<int>::const_iterator p1
=intersection
.begin();
541 QSet
<int>::const_iterator p2
=p1
+1;
542 vtkIdType intersection1
=_nodes
[*p1
];
543 vtkIdType intersection2
=_nodes
[*p2
];
544 if(n2n
[intersection1
].contains(intersection2
))//if there's an edge between intersection1 and intersection2
546 //check if (node1,intersection1,intersection2) and (node2,intersection1,intersection2) are defined as cells!
547 // QVector< QSet< int > > n2c
548 QSet
< int > S1
=n2c
[intersection1
];
549 QSet
< int > S2
=n2c
[intersection2
];
550 QSet
< int > Si
=S1
.intersect(S2
);
552 foreach(vtkIdType C
,Si
){
553 vtkIdType N_pts
, *pts
;
554 grid
->GetCellPoints(C
, N_pts
, pts
);
555 for(int i
=0;i
<N_pts
;i
++)
557 if(pts
[i
]==node1
|| pts
[i
]==node2
) counter
++;
560 if(counter
>=2) IsTetra
=true;
566 // vtkIdType Operation::FindSnapPoint(vtkUnstructuredGrid *src, vtkIdType DeadNode)
571 bool Operation::DeletePoint(vtkUnstructuredGrid
*src
, vtkIdType DeadNode
)
573 EG_VTKSP(vtkUnstructuredGrid
, dst
);
576 int N_points
=src
->GetNumberOfPoints();
577 int N_cells
=src
->GetNumberOfCells();
581 if(DeadNode
<0 || DeadNode
>=N_points
)
583 cout
<<"Warning: Point out of range: DeadNode="<<DeadNode
<<" N_points="<<N_points
<<endl
;
587 QSet
<vtkIdType
> DeadCells
;
588 QSet
<vtkIdType
> MutatedCells
;
589 QSet
<vtkIdType
> MutilatedCells
;
591 vtkIdType SnapPoint
=-1;
592 //Find closest point to DeadNode
593 // vtkIdType SnapPoint = getClosestNode(DeadNode,src);//DeadNode moves to SnapPoint
595 foreach(vtkIdType PSP
, n2n
[DeadNode
])
597 bool IsValidSnapPoint
=true;
599 cout
<<"====>PSP="<<PSP
<<endl
;
601 if(NumberOfCommonPoints(DeadNode
,PSP
,IsTetra
)>2)//common point check
603 cout
<<"Sorry, but you are not allowed to move point "<<DeadNode
<<" to point "<<PSP
<<"."<<endl
;
604 IsValidSnapPoint
=false;
606 if(IsTetra
)//tetra check
608 cout
<<"Sorry, but you are not allowed to move point "<<DeadNode
<<" to point "<<PSP
<<"."<<endl
;
609 IsValidSnapPoint
=false;
612 //count number of points and cells to remove + analyse cell transformations
616 MutatedCells
.clear();
617 MutilatedCells
.clear();
618 foreach(vtkIdType C
, n2c
[DeadNode
])//loop through potentially dead cells
621 //get points around cell
622 vtkIdType N_pts
, *pts
;
623 src
->GetCellPoints(C
, N_pts
, pts
);
625 bool ContainsSnapPoint
=false;
626 bool invincible
=false;
627 for(int i
=0;i
<N_pts
;i
++)
629 cout
<<"pts["<<i
<<"]="<<pts
[i
]<<" and PSP="<<PSP
<<endl
;
630 if(pts
[i
]==PSP
) {ContainsSnapPoint
=true;}
631 if(pts
[i
]!=DeadNode
&& pts
[i
]!=PSP
&& n2c
[pts
[i
]].size()<=1) invincible
=true;
633 if(ContainsSnapPoint
)
635 if(N_pts
==3)//dead cell
637 if(invincible
)//Check that empty lines aren't left behind when a cell is killed
639 cout
<<"Sorry, but you are not allowed to move point "<<DeadNode
<<" to point "<<PSP
<<"."<<endl
;
640 IsValidSnapPoint
=false;
646 cout
<<"cell "<<C
<<" has been pwned!"<<endl
;
649 /* else if(N_pts==4)//mutilated cell
651 MutilatedCells.insert(C);
652 cout<<"cell "<<C<<" has lost a limb!"<<endl;
656 cout
<<"RED ALERT: Xenomorph detected!"<<endl
;
662 vtkIdType src_N_pts
, *src_pts
;
663 src
->GetCellPoints(C
, src_N_pts
, src_pts
);
667 cout
<<"RED ALERT: Xenomorph detected!"<<endl
;
671 vtkIdType OldTriangle
[3];
672 vtkIdType NewTriangle
[3];
674 for(int i
=0;i
<src_N_pts
;i
++)
676 OldTriangle
[i
]=src_pts
[i
];
677 NewTriangle
[i
]=( (src_pts
[i
]==DeadNode
) ? PSP
: src_pts
[i
] );
679 vec3_t Old_N
= triNormal(src
, OldTriangle
[0], OldTriangle
[1], OldTriangle
[2]);
680 vec3_t New_N
= triNormal(src
, NewTriangle
[0], NewTriangle
[1], NewTriangle
[2]);
681 double OldArea
=Old_N
.abs();
682 double NewArea
=New_N
.abs();
683 double scal
=Old_N
*New_N
;
684 double cross
=(Old_N
.cross(New_N
)).abs();//double-cross on Nar Shadaa B-)
686 cout
<<"OldArea="<<OldArea
<<endl
;
687 cout
<<"NewArea="<<NewArea
<<endl
;
688 cout
<<"scal="<<scal
<<endl
;
689 cout
<<"cross="<<cross
<<endl
;
691 if(Old_N
*New_N
<Old_N
*Old_N
*1./100.)//area + inversion check
693 cout
<<"Sorry, but you are not allowed to move point "<<DeadNode
<<" to point "<<PSP
<<"."<<endl
;
694 IsValidSnapPoint
=false;
697 /* if(NewArea<OldArea*1./100.)
699 cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
700 IsValidSnapPoint=false;
705 cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
706 IsValidSnapPoint=false;
710 MutatedCells
.insert(C
);
711 cout
<<"cell "<<C
<<" has been infected!"<<endl
;
715 if(N_cells
+N_newcells
<=0)//survivor check
717 cout
<<"Sorry, but you are not allowed to move point "<<DeadNode
<<" to point "<<PSP
<<"."<<endl
;
718 IsValidSnapPoint
=false;
720 /* if(EmptyVolume(DeadNode,PSP))//simplified volume check
722 cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
723 IsValidSnapPoint=false;
725 if(IsValidSnapPoint
) {SnapPoint
=PSP
; break;}
726 }//end of loop through potential SnapPoints
728 cout
<<"===>SNAPPOINT="<<SnapPoint
<<endl
;
729 if(SnapPoint
<0) {cout
<<"Sorry no possible SnapPoint found."<<endl
; return(false);}
732 cout
<<"N_points="<<N_points
<<endl
;
733 cout
<<"N_cells="<<N_cells
<<endl
;
734 cout
<<"N_newpoints="<<N_newpoints
<<endl
;
735 cout
<<"N_newcells="<<N_newcells
<<endl
;
736 allocateGrid(dst
,N_cells
+N_newcells
,N_points
+N_newpoints
);
738 //vector used to redefine the new point IDs
739 QVector
<vtkIdType
> OffSet(N_points
);
742 vtkIdType dst_id_node
=0;
743 for (vtkIdType src_id_node
= 0; src_id_node
< N_points
; src_id_node
++) {//loop through src points
744 if(src_id_node
!=DeadNode
)//if the node isn't dead, copy it
747 src
->GetPoints()->GetPoint(src_id_node
, x
.data());
748 dst
->GetPoints()->SetPoint(dst_id_node
, x
.data());
749 copyNodeData(src
, src_id_node
, dst
, dst_id_node
);
750 OffSet
[src_id_node
]=src_id_node
-dst_id_node
;
755 cout
<<"DeadCells="<<DeadCells
<<endl
;
756 cout
<<"MutatedCells="<<MutatedCells
<<endl
;
757 cout
<<"MutilatedCells="<<MutilatedCells
<<endl
;
760 for (vtkIdType id_cell
= 0; id_cell
< src
->GetNumberOfCells(); ++id_cell
) {//loop through src cells
761 if(!DeadCells
.contains(id_cell
))//if the cell isn't dead
763 vtkIdType src_N_pts
, *src_pts
;
764 vtkIdType dst_N_pts
, *dst_pts
;
765 src
->GetCellPoints(id_cell
, src_N_pts
, src_pts
);
767 vtkIdType type_cell
= src
->GetCellType(id_cell
);
768 cout
<<"-->id_cell="<<id_cell
<<endl
;
769 for(int i
=0;i
<src_N_pts
;i
++) cout
<<"src_pts["<<i
<<"]="<<src_pts
[i
]<<endl
;
770 // src->GetCellPoints(id_cell, dst_N_pts, dst_pts);
772 dst_pts
=new vtkIdType
[dst_N_pts
];
773 if(MutatedCells
.contains(id_cell
))//mutated cell
775 cout
<<"processing mutated cell "<<id_cell
<<endl
;
776 for(int i
=0;i
<src_N_pts
;i
++)
778 if(src_pts
[i
]==DeadNode
) {
779 cout
<<"SnapPoint="<<SnapPoint
<<endl
;
780 cout
<<"OffSet[SnapPoint]="<<OffSet
[SnapPoint
]<<endl
;
781 cout
<<"src_pts["<<i
<<"]="<<src_pts
[i
]<<endl
;
782 dst_pts
[i
]=SnapPoint
-OffSet
[SnapPoint
];
784 else dst_pts
[i
]=src_pts
[i
]-OffSet
[src_pts
[i
]];
786 cout
<<"--->dst_pts:"<<endl
;
787 for(int i
=0;i
<dst_N_pts
;i
++) cout
<<"dst_pts["<<i
<<"]="<<dst_pts
[i
]<<endl
;
790 else if(MutilatedCells
.contains(id_cell
))//mutilated cell
792 cout
<<"processing mutilated cell "<<id_cell
<<endl
;
794 if(type_cell
==VTK_QUAD
) {
795 type_cell
=VTK_TRIANGLE
;
798 else {cout
<<"FATAL ERROR: Unknown mutilated cell detected! It is not a quad! Potential xenomorph infestation!"<<endl
;EG_BUG
;}
801 for(int i
=0;i
<src_N_pts
;i
++)
803 if(src_pts
[i
]==SnapPoint
) { dst_pts
[j
]=SnapPoint
-OffSet
[SnapPoint
];j
++; }//SnapPoint
804 else if(src_pts
[i
]!=DeadNode
) { dst_pts
[j
]=src_pts
[i
]-OffSet
[src_pts
[i
]];j
++; }//pre-snap/dead + post-snap/dead
805 //do nothing in case of DeadNode
810 cout
<<"processing normal cell "<<id_cell
<<endl
;
811 for(int i
=0;i
<src_N_pts
;i
++)
813 dst_pts
[i
]=src_pts
[i
]-OffSet
[src_pts
[i
]];
817 vtkIdType id_new_cell
= dst
->InsertNextCell(type_cell
, dst_N_pts
, dst_pts
);
818 copyCellData(src
, id_cell
, dst
, id_new_cell
);
819 cout
<<"===Copying cell "<<id_cell
<<" to "<<id_new_cell
<<"==="<<endl
;
820 cout
<<"src_pts:"<<endl
;
821 for(int i
=0;i
<src_N_pts
;i
++) cout
<<"src_pts["<<i
<<"]="<<src_pts
[i
]<<endl
;
822 cout
<<"dst_pts:"<<endl
;
823 for(int i
=0;i
<dst_N_pts
;i
++) cout
<<"dst_pts["<<i
<<"]="<<dst_pts
[i
]<<endl
;
824 cout
<<"OffSet="<<OffSet
<<endl
;
825 cout
<<"===Copying cell end==="<<endl
;
829 // cout_grid(cout,dst,true,true,true,true);
834 bool Operation::EmptyVolume(vtkIdType DeadNode
, vtkIdType PSP
)
841 vec3_t
Operation::GetCenter(vtkIdType cellId
, double& R
)
843 vtkIdType
*pts
, Npts
;
844 grid
->GetCellPoints(cellId
, Npts
, pts
);
847 for (vtkIdType i
= 0; i
< Npts
; ++i
) {
849 grid
->GetPoints()->GetPoint(pts
[i
], xp
.data());
850 x
+= double(1)/Npts
* xp
;
854 for (vtkIdType i
= 0; i
< Npts
; ++i
) {
856 grid
->GetPoints()->GetPoint(pts
[i
], xp
.data());
857 R
= min(R
, 0.25*(xp
-x
).abs());