1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkMolecule.cxx
5 Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
6 All rights reserved.
7 See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
8
9 This software is distributed WITHOUT ANY WARRANTY; without even
10 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
11 PURPOSE. See the above copyright notice for more information.
12
13 =========================================================================*/
14 #include "vtkMolecule.h"
15
16 #include "vtkAbstractElectronicData.h"
17 #include "vtkDataSetAttributes.h"
18 #include "vtkEdgeListIterator.h"
19 #include "vtkFloatArray.h"
20 #include "vtkGraphInternals.h"
21 #include "vtkIdTypeArray.h"
22 #include "vtkInformation.h"
23 #include "vtkInformationVector.h"
24 #include "vtkMatrix3x3.h"
25 #include "vtkNew.h"
26 #include "vtkObjectFactory.h"
27 #include "vtkPlane.h"
28 #include "vtkPoints.h"
29 #include "vtkUnsignedCharArray.h"
30 #include "vtkUnsignedShortArray.h"
31 #include "vtkVector.h"
32 #include "vtkVectorOperators.h"
33
34 #include <cassert>
35
36 //------------------------------------------------------------------------------
37 vtkStandardNewMacro(vtkMolecule);
38
39 //------------------------------------------------------------------------------
vtkMolecule()40 vtkMolecule::vtkMolecule()
41 : ElectronicData(nullptr)
42 , Lattice(nullptr)
43 , LatticeOrigin(0., 0., 0.)
44 , AtomGhostArray(nullptr)
45 , BondGhostArray(nullptr)
46 , AtomicNumberArrayName(nullptr)
47 , BondOrdersArrayName(nullptr)
48 {
49 this->Initialize();
50 }
51
52 //------------------------------------------------------------------------------
Initialize()53 void vtkMolecule::Initialize()
54 {
55 // Reset underlying data structure
56 this->Superclass::Initialize();
57
58 // Setup vertex data
59 vtkDataSetAttributes* vertData = this->GetVertexData();
60 vertData->AllocateArrays(1); // atomic nums
61
62 // Atomic numbers
63 this->SetAtomicNumberArrayName("Atomic Numbers");
64 vtkNew<vtkUnsignedShortArray> atomicNums;
65 atomicNums->SetNumberOfComponents(1);
66 atomicNums->SetName(this->GetAtomicNumberArrayName());
67 vertData->SetScalars(atomicNums);
68
69 // Nuclear coordinates
70 vtkPoints* points = vtkPoints::New();
71 this->SetPoints(points);
72 points->Delete();
73
74 // Setup edge data
75 vtkDataSetAttributes* edgeData = this->GetEdgeData();
76 edgeData->AllocateArrays(1); // Bond orders
77
78 this->SetBondOrdersArrayName("Bond Orders");
79 vtkNew<vtkUnsignedShortArray> bondOrders;
80 bondOrders->SetNumberOfComponents(1);
81 bondOrders->SetName(this->GetBondOrdersArrayName());
82 edgeData->SetScalars(bondOrders);
83
84 this->UpdateBondList();
85
86 // Electronic data
87 this->SetElectronicData(nullptr);
88
89 this->Modified();
90 }
91
92 //------------------------------------------------------------------------------
~vtkMolecule()93 vtkMolecule::~vtkMolecule()
94 {
95 this->SetElectronicData(nullptr);
96 delete[] this->AtomicNumberArrayName;
97 delete[] this->BondOrdersArrayName;
98 }
99
100 //------------------------------------------------------------------------------
PrintSelf(ostream & os,vtkIndent indent)101 void vtkMolecule::PrintSelf(ostream& os, vtkIndent indent)
102 {
103 this->Superclass::PrintSelf(os, indent);
104
105 vtkIndent subIndent = indent.GetNextIndent();
106
107 os << indent << "Atoms:\n";
108 for (vtkIdType i = 0; i < this->GetNumberOfAtoms(); ++i)
109 {
110 this->GetAtom(i).PrintSelf(os, subIndent);
111 }
112
113 os << indent << "Bonds:\n";
114 for (vtkIdType i = 0; i < this->GetNumberOfBonds(); ++i)
115 {
116 os << subIndent << "===== Bond " << i << ": =====\n";
117 this->GetBond(i).PrintSelf(os, subIndent);
118 }
119
120 os << indent << "Lattice:\n";
121 if (this->HasLattice())
122 {
123 double* m = this->Lattice->GetData();
124 os << subIndent << "a: " << m[0] << " " << m[3] << " " << m[6] << "\n";
125 os << subIndent << "b: " << m[1] << " " << m[4] << " " << m[7] << "\n";
126 os << subIndent << "c: " << m[2] << " " << m[5] << " " << m[8] << "\n";
127 os << subIndent << "origin: " << this->LatticeOrigin[0] << " " << this->LatticeOrigin[1] << " "
128 << this->LatticeOrigin[2] << "\n";
129 }
130
131 os << indent << "Electronic Data:\n";
132 if (this->ElectronicData)
133 {
134 this->ElectronicData->PrintSelf(os, subIndent);
135 }
136 else
137 {
138 os << subIndent << "Not set.\n";
139 }
140
141 os << indent << "Atomic number array name : " << this->GetAtomicNumberArrayName() << "\n";
142 os << indent << "Bond orders array name : " << this->GetBondOrdersArrayName();
143 }
144
145 //------------------------------------------------------------------------------
AppendAtom(unsigned short atomicNumber,double x,double y,double z)146 vtkAtom vtkMolecule::AppendAtom(unsigned short atomicNumber, double x, double y, double z)
147 {
148 vtkUnsignedShortArray* atomicNums = this->GetAtomicNumberArray();
149
150 assert(atomicNums);
151
152 vtkIdType id;
153 this->AddVertexInternal(nullptr, &id);
154
155 atomicNums->InsertValue(id, atomicNumber);
156 vtkIdType coordID = this->Points->InsertNextPoint(x, y, z);
157 (void)coordID;
158 assert("point ids synced with vertex ids" && coordID == id);
159
160 this->Modified();
161 return vtkAtom(this, id);
162 }
163
164 //------------------------------------------------------------------------------
GetAtom(vtkIdType atomId)165 vtkAtom vtkMolecule::GetAtom(vtkIdType atomId)
166 {
167 assert(atomId >= 0 && atomId < this->GetNumberOfAtoms());
168
169 vtkAtom atom(this, atomId);
170 return atom;
171 }
172
173 //------------------------------------------------------------------------------
GetAtomAtomicNumber(vtkIdType id)174 unsigned short vtkMolecule::GetAtomAtomicNumber(vtkIdType id)
175 {
176 assert(id >= 0 && id < this->GetNumberOfAtoms());
177
178 vtkUnsignedShortArray* atomicNums = this->GetAtomicNumberArray();
179
180 return atomicNums->GetValue(id);
181 }
182
183 //------------------------------------------------------------------------------
SetAtomAtomicNumber(vtkIdType id,unsigned short atomicNum)184 void vtkMolecule::SetAtomAtomicNumber(vtkIdType id, unsigned short atomicNum)
185 {
186 assert(id >= 0 && id < this->GetNumberOfAtoms());
187
188 vtkUnsignedShortArray* atomicNums = this->GetAtomicNumberArray();
189
190 atomicNums->SetValue(id, atomicNum);
191 this->Modified();
192 }
193
194 //------------------------------------------------------------------------------
SetAtomPosition(vtkIdType id,const vtkVector3f & pos)195 void vtkMolecule::SetAtomPosition(vtkIdType id, const vtkVector3f& pos)
196 {
197 assert(id >= 0 && id < this->GetNumberOfAtoms());
198 this->Points->SetPoint(id, pos.GetData());
199 this->Modified();
200 }
201
202 //------------------------------------------------------------------------------
SetAtomPosition(vtkIdType id,double x,double y,double z)203 void vtkMolecule::SetAtomPosition(vtkIdType id, double x, double y, double z)
204 {
205 assert(id >= 0 && id < this->GetNumberOfAtoms());
206 this->Points->SetPoint(id, x, y, z);
207 this->Modified();
208 }
209
210 //------------------------------------------------------------------------------
GetAtomPosition(vtkIdType id)211 vtkVector3f vtkMolecule::GetAtomPosition(vtkIdType id)
212 {
213 assert(id >= 0 && id < this->GetNumberOfAtoms());
214 vtkDataArray* positions = this->Points->GetData();
215 auto positionsF = vtkArrayDownCast<vtkFloatArray>(positions);
216 if (positionsF)
217 {
218 float* data = positionsF->GetPointer(id * 3);
219 return vtkVector3f(data);
220 }
221
222 auto point = positions->GetTuple3(id);
223 return vtkVector3f(point[0], point[1], point[2]);
224 }
225
226 //------------------------------------------------------------------------------
GetAtomPosition(vtkIdType id,float pos[3])227 void vtkMolecule::GetAtomPosition(vtkIdType id, float pos[3])
228 {
229 vtkVector3f position = this->GetAtomPosition(id);
230 pos[0] = position.GetX();
231 pos[1] = position.GetY();
232 pos[2] = position.GetZ();
233 }
234
235 //------------------------------------------------------------------------------
GetAtomPosition(vtkIdType id,double pos[3])236 void vtkMolecule::GetAtomPosition(vtkIdType id, double pos[3])
237 {
238 this->Points->GetPoint(id, pos);
239 }
240
241 //------------------------------------------------------------------------------
GetNumberOfAtoms()242 vtkIdType vtkMolecule::GetNumberOfAtoms()
243 {
244 return this->GetNumberOfVertices();
245 }
246
247 //------------------------------------------------------------------------------
AppendBond(const vtkIdType atom1,const vtkIdType atom2,const unsigned short order)248 vtkBond vtkMolecule::AppendBond(
249 const vtkIdType atom1, const vtkIdType atom2, const unsigned short order)
250 {
251 vtkUnsignedShortArray* bondOrders = this->GetBondOrdersArray();
252
253 assert(bondOrders);
254
255 vtkEdgeType edgeType;
256 this->AddEdgeInternal(atom1, atom2, false, nullptr, &edgeType);
257 this->SetBondListDirty();
258
259 vtkIdType id = edgeType.Id;
260 bondOrders->InsertValue(id, order);
261 this->Modified();
262 return vtkBond(this, id, atom1, atom2);
263 }
264
265 //------------------------------------------------------------------------------
GetBond(vtkIdType bondId)266 vtkBond vtkMolecule::GetBond(vtkIdType bondId)
267 {
268 assert(bondId >= 0 && bondId < this->GetNumberOfBonds());
269
270 vtkIdTypeArray* bonds = this->GetBondList();
271 // An array with two components holding the bonded atom's ids
272 vtkIdType* ids = bonds->GetPointer(2 * bondId);
273 return vtkBond(this, bondId, ids[0], ids[1]);
274 }
275
276 //------------------------------------------------------------------------------
SetBondOrder(vtkIdType bondId,unsigned short order)277 void vtkMolecule::SetBondOrder(vtkIdType bondId, unsigned short order)
278 {
279 assert(bondId >= 0 && bondId < this->GetNumberOfBonds());
280
281 vtkUnsignedShortArray* bondOrders = this->GetBondOrdersArray();
282 assert(bondOrders);
283
284 this->Modified();
285 bondOrders->InsertValue(bondId, order);
286 }
287
288 //------------------------------------------------------------------------------
GetBondOrder(vtkIdType bondId)289 unsigned short vtkMolecule::GetBondOrder(vtkIdType bondId)
290 {
291 assert(bondId >= 0 && bondId < this->GetNumberOfBonds());
292
293 vtkUnsignedShortArray* bondOrders = this->GetBondOrdersArray();
294
295 return bondOrders ? bondOrders->GetValue(bondId) : 0;
296 }
297
298 //------------------------------------------------------------------------------
GetBondLength(vtkIdType bondId)299 double vtkMolecule::GetBondLength(vtkIdType bondId)
300 {
301 vtkBond bond = this->GetBond(bondId);
302 return bond.GetLength();
303 }
304
305 //------------------------------------------------------------------------------
GetAtomicPositionArray()306 vtkPoints* vtkMolecule::GetAtomicPositionArray()
307 {
308 return this->Points;
309 }
310
311 //------------------------------------------------------------------------------
GetAtomicNumberArray()312 vtkUnsignedShortArray* vtkMolecule::GetAtomicNumberArray()
313 {
314 vtkUnsignedShortArray* atomicNums = vtkArrayDownCast<vtkUnsignedShortArray>(
315 this->GetVertexData()->GetScalars(this->GetAtomicNumberArrayName()));
316
317 assert(atomicNums);
318
319 return atomicNums;
320 }
321
322 //------------------------------------------------------------------------------
GetBondOrdersArray()323 vtkUnsignedShortArray* vtkMolecule::GetBondOrdersArray()
324 {
325 return vtkArrayDownCast<vtkUnsignedShortArray>(
326 this->GetBondData()->GetScalars(this->GetBondOrdersArrayName()));
327 }
328
329 //------------------------------------------------------------------------------
GetNumberOfBonds()330 vtkIdType vtkMolecule::GetNumberOfBonds()
331 {
332 return this->GetNumberOfEdges();
333 }
334
335 //------------------------------------------------------------------------------
336 vtkCxxSetObjectMacro(vtkMolecule, ElectronicData, vtkAbstractElectronicData);
337
338 //------------------------------------------------------------------------------
ShallowCopy(vtkDataObject * obj)339 void vtkMolecule::ShallowCopy(vtkDataObject* obj)
340 {
341 vtkMolecule* m = vtkMolecule::SafeDownCast(obj);
342 if (!m)
343 {
344 vtkErrorMacro("Can only shallow copy from vtkMolecule or subclass.");
345 return;
346 }
347 this->ShallowCopyStructure(m);
348 this->ShallowCopyAttributes(m);
349 }
350
351 //------------------------------------------------------------------------------
DeepCopy(vtkDataObject * obj)352 void vtkMolecule::DeepCopy(vtkDataObject* obj)
353 {
354 vtkMolecule* m = vtkMolecule::SafeDownCast(obj);
355 if (!m)
356 {
357 vtkErrorMacro("Can only deep copy from vtkMolecule or subclass.");
358 return;
359 }
360 this->DeepCopyStructure(m);
361 this->DeepCopyAttributes(m);
362 }
363
364 //------------------------------------------------------------------------------
CheckedShallowCopy(vtkGraph * g)365 bool vtkMolecule::CheckedShallowCopy(vtkGraph* g)
366 {
367 bool result = this->Superclass::CheckedShallowCopy(g);
368 this->BondListIsDirty = true;
369 return result;
370 }
371
372 //------------------------------------------------------------------------------
CheckedDeepCopy(vtkGraph * g)373 bool vtkMolecule::CheckedDeepCopy(vtkGraph* g)
374 {
375 bool result = this->Superclass::CheckedDeepCopy(g);
376 this->BondListIsDirty = true;
377 return result;
378 }
379
380 //------------------------------------------------------------------------------
ShallowCopyStructure(vtkMolecule * m)381 void vtkMolecule::ShallowCopyStructure(vtkMolecule* m)
382 {
383 this->CopyStructureInternal(m, false);
384 }
385
386 //------------------------------------------------------------------------------
DeepCopyStructure(vtkMolecule * m)387 void vtkMolecule::DeepCopyStructure(vtkMolecule* m)
388 {
389 this->CopyStructureInternal(m, true);
390 }
391
392 //------------------------------------------------------------------------------
ShallowCopyAttributes(vtkMolecule * m)393 void vtkMolecule::ShallowCopyAttributes(vtkMolecule* m)
394 {
395 this->CopyAttributesInternal(m, false);
396 }
397
398 //------------------------------------------------------------------------------
DeepCopyAttributes(vtkMolecule * m)399 void vtkMolecule::DeepCopyAttributes(vtkMolecule* m)
400 {
401 this->CopyAttributesInternal(m, true);
402 }
403
404 //------------------------------------------------------------------------------
CopyStructureInternal(vtkMolecule * m,bool deep)405 void vtkMolecule::CopyStructureInternal(vtkMolecule* m, bool deep)
406 {
407 // Call superclass
408 if (deep)
409 {
410 this->Superclass::DeepCopy(m);
411 }
412 else
413 {
414 this->Superclass::ShallowCopy(m);
415 }
416
417 if (!m->HasLattice())
418 {
419 this->ClearLattice();
420 }
421 else
422 {
423 if (deep)
424 {
425 vtkNew<vtkMatrix3x3> newLattice;
426 newLattice->DeepCopy(m->Lattice);
427 this->SetLattice(newLattice);
428 }
429 else
430 {
431 this->SetLattice(m->Lattice);
432 }
433 this->LatticeOrigin = m->LatticeOrigin;
434 }
435
436 this->BondListIsDirty = true;
437 }
438
439 //------------------------------------------------------------------------------
CopyAttributesInternal(vtkMolecule * m,bool deep)440 void vtkMolecule::CopyAttributesInternal(vtkMolecule* m, bool deep)
441 {
442 if (deep)
443 {
444 if (m->ElectronicData)
445 this->ElectronicData->DeepCopy(m->ElectronicData);
446 }
447 else
448 {
449 this->SetElectronicData(m->ElectronicData);
450 }
451 }
452
453 //------------------------------------------------------------------------------
UpdateBondList()454 void vtkMolecule::UpdateBondList()
455 {
456 this->BuildEdgeList();
457 this->BondListIsDirty = false;
458 }
459
460 //------------------------------------------------------------------------------
GetBondList()461 vtkIdTypeArray* vtkMolecule::GetBondList()
462 {
463 // Create the edge list if it doesn't exist, or is marked as dirty.
464 vtkIdTypeArray* edgeList = this->BondListIsDirty ? nullptr : this->GetEdgeList();
465 if (!edgeList)
466 {
467 this->UpdateBondList();
468 edgeList = this->GetEdgeList();
469 }
470 assert(edgeList != nullptr);
471 return edgeList;
472 }
473
474 //------------------------------------------------------------------------------
GetPlaneFromBond(const vtkBond & bond,const vtkVector3f & normal,vtkPlane * plane)475 bool vtkMolecule::GetPlaneFromBond(const vtkBond& bond, const vtkVector3f& normal, vtkPlane* plane)
476 {
477 return vtkMolecule::GetPlaneFromBond(bond.GetBeginAtom(), bond.GetEndAtom(), normal, plane);
478 }
479
480 //------------------------------------------------------------------------------
GetPlaneFromBond(const vtkAtom & atom1,const vtkAtom & atom2,const vtkVector3f & normal,vtkPlane * plane)481 bool vtkMolecule::GetPlaneFromBond(
482 const vtkAtom& atom1, const vtkAtom& atom2, const vtkVector3f& normal, vtkPlane* plane)
483 {
484 if (plane == nullptr)
485 {
486 return false;
487 }
488
489 vtkVector3f v(atom1.GetPosition() - atom2.GetPosition());
490
491 vtkVector3f n_i(normal);
492 vtkVector3f unitV(v.Normalized());
493
494 // Check if vectors are (nearly) parallel
495 if (unitV.Compare(n_i.Normalized(), 1e-7))
496 {
497 return false;
498 }
499
500 // calculate projection of n_i onto v
501 // TODO Remove or restore this when scalar mult. is supported again
502 // vtkVector3d proj (unitV * n_i.Dot(unitV));
503 double n_iDotUnitV = n_i.Dot(unitV);
504 vtkVector3f proj(unitV[0] * n_iDotUnitV, unitV[1] * n_iDotUnitV, unitV[2] * n_iDotUnitV);
505 // end vtkVector reimplementation TODO
506
507 // Calculate actual normal:
508 vtkVector3f realNormal(n_i - proj);
509
510 // Create plane:
511 vtkVector3f pos(atom1.GetPosition());
512 plane->SetOrigin(pos.Cast<double>().GetData());
513 plane->SetNormal(realNormal.Cast<double>().GetData());
514 return true;
515 }
516
517 //------------------------------------------------------------------------------
HasLattice()518 bool vtkMolecule::HasLattice()
519 {
520 return this->Lattice != nullptr;
521 }
522
523 //------------------------------------------------------------------------------
ClearLattice()524 void vtkMolecule::ClearLattice()
525 {
526 this->SetLattice(nullptr);
527 }
528
529 //------------------------------------------------------------------------------
SetLattice(vtkMatrix3x3 * matrix)530 void vtkMolecule::SetLattice(vtkMatrix3x3* matrix)
531 {
532 if (!matrix)
533 {
534 if (this->Lattice)
535 {
536 // If we're clearing a matrix, zero out the origin:
537 this->LatticeOrigin = vtkVector3d(0., 0., 0.);
538 this->Lattice = nullptr;
539 this->Modified();
540 }
541 }
542 else if (this->Lattice != matrix)
543 {
544 this->Lattice = matrix;
545 this->Modified();
546 }
547 }
548
549 //------------------------------------------------------------------------------
SetLattice(const vtkVector3d & a,const vtkVector3d & b,const vtkVector3d & c)550 void vtkMolecule::SetLattice(const vtkVector3d& a, const vtkVector3d& b, const vtkVector3d& c)
551 {
552 if (this->Lattice == nullptr)
553 {
554 this->Lattice.TakeReference(vtkMatrix3x3::New());
555 this->Modified();
556 }
557
558 double* mat = this->Lattice->GetData();
559 if (mat[0] != a[0] || mat[1] != b[0] || mat[2] != c[0] || mat[3] != a[1] || mat[4] != b[1] ||
560 mat[5] != c[1] || mat[6] != a[2] || mat[7] != b[2] || mat[8] != c[2])
561 {
562 mat[0] = a[0];
563 mat[1] = b[0];
564 mat[2] = c[0];
565 mat[3] = a[1];
566 mat[4] = b[1];
567 mat[5] = c[1];
568 mat[6] = a[2];
569 mat[7] = b[2];
570 mat[8] = c[2];
571 this->Modified();
572 }
573 }
574
575 //------------------------------------------------------------------------------
GetLattice()576 vtkMatrix3x3* vtkMolecule::GetLattice()
577 {
578 return this->Lattice;
579 }
580
581 //------------------------------------------------------------------------------
GetLattice(vtkVector3d & a,vtkVector3d & b,vtkVector3d & c)582 void vtkMolecule::GetLattice(vtkVector3d& a, vtkVector3d& b, vtkVector3d& c)
583 {
584 if (this->Lattice)
585 {
586 double* mat = this->Lattice->GetData();
587 a[0] = mat[0];
588 a[1] = mat[3];
589 a[2] = mat[6];
590 b[0] = mat[1];
591 b[1] = mat[4];
592 b[2] = mat[7];
593 c[0] = mat[2];
594 c[1] = mat[5];
595 c[2] = mat[8];
596 }
597 else
598 {
599 a = b = c = vtkVector3d(0., 0., 0.);
600 }
601 }
602
603 //------------------------------------------------------------------------------
GetLattice(vtkVector3d & a,vtkVector3d & b,vtkVector3d & c,vtkVector3d & origin)604 void vtkMolecule::GetLattice(vtkVector3d& a, vtkVector3d& b, vtkVector3d& c, vtkVector3d& origin)
605 {
606 if (this->Lattice)
607 {
608 double* mat = this->Lattice->GetData();
609 a[0] = mat[0];
610 a[1] = mat[3];
611 a[2] = mat[6];
612 b[0] = mat[1];
613 b[1] = mat[4];
614 b[2] = mat[7];
615 c[0] = mat[2];
616 c[1] = mat[5];
617 c[2] = mat[8];
618 origin = this->LatticeOrigin;
619 }
620 else
621 {
622 a = b = c = origin = vtkVector3d(0., 0., 0.);
623 }
624 }
625
626 //------------------------------------------------------------------------------
GetAtomGhostArray()627 vtkUnsignedCharArray* vtkMolecule::GetAtomGhostArray()
628 {
629 return vtkArrayDownCast<vtkUnsignedCharArray>(
630 this->GetVertexData()->GetArray(vtkDataSetAttributes::GhostArrayName()));
631 }
632
633 //------------------------------------------------------------------------------
AllocateAtomGhostArray()634 void vtkMolecule::AllocateAtomGhostArray()
635 {
636 if (this->GetAtomGhostArray() == nullptr)
637 {
638 vtkNew<vtkUnsignedCharArray> ghosts;
639 ghosts->SetName(vtkDataSetAttributes::GhostArrayName());
640 ghosts->SetNumberOfComponents(1);
641 ghosts->SetNumberOfTuples(this->GetNumberOfAtoms());
642 ghosts->FillComponent(0, 0);
643 this->GetVertexData()->AddArray(ghosts);
644 }
645 else
646 {
647 this->GetAtomGhostArray()->SetNumberOfTuples(this->GetNumberOfAtoms());
648 }
649 }
650
651 //------------------------------------------------------------------------------
GetBondGhostArray()652 vtkUnsignedCharArray* vtkMolecule::GetBondGhostArray()
653 {
654 return vtkArrayDownCast<vtkUnsignedCharArray>(
655 this->GetEdgeData()->GetArray(vtkDataSetAttributes::GhostArrayName()));
656 }
657
658 //------------------------------------------------------------------------------
AllocateBondGhostArray()659 void vtkMolecule::AllocateBondGhostArray()
660 {
661 if (this->GetBondGhostArray() == nullptr)
662 {
663 vtkNew<vtkUnsignedCharArray> ghosts;
664 ghosts->SetName(vtkDataSetAttributes::GhostArrayName());
665 ghosts->SetNumberOfComponents(1);
666 ghosts->SetNumberOfTuples(this->GetNumberOfBonds());
667 ghosts->FillComponent(0, 0);
668 this->GetEdgeData()->AddArray(ghosts);
669 }
670 else
671 {
672 this->GetBondGhostArray()->SetNumberOfTuples(this->GetNumberOfBonds());
673 }
674 }
675
676 //------------------------------------------------------------------------------
Initialize(vtkPoints * atomPositions,vtkDataArray * atomicNumberArray,vtkDataSetAttributes * atomData)677 int vtkMolecule::Initialize(
678 vtkPoints* atomPositions, vtkDataArray* atomicNumberArray, vtkDataSetAttributes* atomData)
679 {
680 // Start with default initialization the molecule
681 this->Initialize();
682
683 // if no atomicNumberArray given, look for one in atomData
684 if (!atomicNumberArray && atomData)
685 {
686 atomicNumberArray = atomData->GetArray(this->GetAtomicNumberArrayName());
687 }
688
689 if (!atomPositions && !atomicNumberArray)
690 {
691 vtkDebugMacro(<< "Atom position and atomic numbers were not found: skip atomic data.");
692 return 1;
693 }
694
695 if (!atomPositions || !atomicNumberArray)
696 {
697 vtkDebugMacro(<< "Empty atoms or atomic numbers.");
698 return 0;
699 }
700
701 // ensure it is a short array
702 vtkNew<vtkUnsignedShortArray> newAtomicNumberShortArray;
703 if (vtkUnsignedShortArray::SafeDownCast(atomicNumberArray))
704 {
705 newAtomicNumberShortArray->ShallowCopy(atomicNumberArray);
706 }
707 else
708 {
709 vtkIdType nbPoints = atomicNumberArray->GetNumberOfTuples();
710 newAtomicNumberShortArray->SetNumberOfComponents(1);
711 newAtomicNumberShortArray->SetNumberOfTuples(nbPoints);
712 newAtomicNumberShortArray->SetName(atomicNumberArray->GetName());
713 for (vtkIdType i = 0; i < nbPoints; i++)
714 {
715 newAtomicNumberShortArray->SetTuple1(i, atomicNumberArray->GetTuple1(i));
716 }
717 }
718
719 int nbAtoms = atomPositions->GetNumberOfPoints();
720 if (nbAtoms != newAtomicNumberShortArray->GetNumberOfTuples())
721 {
722 vtkErrorMacro(<< "Number of atoms not equal to number of atomic numbers.");
723 return 0;
724 }
725 if (atomData && nbAtoms != atomData->GetNumberOfTuples())
726 {
727 vtkErrorMacro(<< "Number of atoms not equal to number of atom properties.");
728 return 0;
729 }
730
731 static const std::string atomicNumberName = this->GetAtomicNumberArrayName();
732
733 // update atoms positions
734 this->ForceOwnership();
735 this->Internals->Adjacency.resize(nbAtoms, vtkVertexAdjacencyList());
736 this->SetPoints(atomPositions);
737
738 // if atom properties exists, copy it in VertexData and look for an atomic number in its arrays.
739 if (atomData)
740 {
741 this->GetVertexData()->ShallowCopy(atomData);
742
743 // if atomData contains an array with the atomic number name,
744 // copy this array with a new name as we will overwrite it.
745 vtkDataArray* otherArray = atomData->GetArray(atomicNumberName.c_str());
746 if (otherArray && otherArray != static_cast<vtkAbstractArray*>(atomicNumberArray))
747 {
748 this->GetVertexData()->RemoveArray(atomicNumberName.c_str());
749
750 // create a new name for the copied array.
751 std::string newName = std::string("Original ") + atomicNumberName;
752
753 // if the new name is available create a copy of the array with another name, and add it
754 // else no backup is done, array will be replaced.
755 if (!atomData->GetArray(newName.c_str()))
756 {
757 vtkDataArray* otherArrayCopy = otherArray->NewInstance();
758 otherArrayCopy->ShallowCopy(otherArray);
759 otherArrayCopy->SetName(newName.c_str());
760 this->GetVertexData()->AddArray(otherArrayCopy);
761 otherArrayCopy->Delete();
762 }
763 else
764 {
765 vtkWarningMacro(<< "Array '" << atomicNumberName << "' was replaced.");
766 }
767 }
768 }
769
770 // add atomic number array: if given array has the expected name, add it directly
771 // (it will replace the current one). Else copy it and add it with atomic number name.
772 if (atomicNumberName == newAtomicNumberShortArray->GetName())
773 {
774 this->GetVertexData()->AddArray(newAtomicNumberShortArray);
775 }
776 else
777 {
778 vtkNew<vtkUnsignedShortArray> atomicNumberArrayCopy;
779 atomicNumberArrayCopy->ShallowCopy(newAtomicNumberShortArray);
780 atomicNumberArrayCopy->SetName(atomicNumberName.c_str());
781 this->GetVertexData()->AddArray(atomicNumberArrayCopy.Get());
782 }
783
784 this->Modified();
785 return 1;
786 }
787
788 //------------------------------------------------------------------------------
Initialize(vtkMolecule * molecule)789 int vtkMolecule::Initialize(vtkMolecule* molecule)
790 {
791 if (molecule == nullptr)
792 {
793 this->Initialize();
794 return 1;
795 }
796
797 return this->Initialize(
798 molecule->GetPoints(), molecule->GetAtomicNumberArray(), molecule->GetVertexData());
799 }
800
801 //------------------------------------------------------------------------------
GetData(vtkInformation * info)802 vtkMolecule* vtkMolecule::GetData(vtkInformation* info)
803 {
804 return info ? vtkMolecule::SafeDownCast(info->Get(DATA_OBJECT())) : nullptr;
805 }
806
807 //------------------------------------------------------------------------------
GetData(vtkInformationVector * v,int i)808 vtkMolecule* vtkMolecule::GetData(vtkInformationVector* v, int i)
809 {
810 return vtkMolecule::GetData(v->GetInformationObject(i));
811 }
812
813 //------------------------------------------------------------------------------
GetActualMemorySize()814 unsigned long vtkMolecule::GetActualMemorySize()
815 {
816 unsigned long size = this->Superclass::GetActualMemorySize();
817 if (this->ElectronicData)
818 {
819 size += this->ElectronicData->GetActualMemorySize();
820 }
821 if (this->AtomGhostArray)
822 {
823 size += this->AtomGhostArray->GetActualMemorySize();
824 }
825 if (this->BondGhostArray)
826 {
827 size += this->BondGhostArray->GetActualMemorySize();
828 }
829 return size;
830 }
831