Changeset a567c3
- Timestamp:
- Jul 28, 2009, 10:33:53 AM (16 years ago)
- Branches:
- Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable
- Children:
- d4d0dd
- Parents:
- fcc9b0
- git-author:
- Frederik Heber <heber@…> (07/28/09 10:31:26)
- git-committer:
- Frederik Heber <heber@…> (07/28/09 10:33:53)
- Location:
- src
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
src/boundary.cpp
rfcc9b0 ra567c3 682 682 void write_tecplot_file(ofstream *out, ofstream *tecplot, class Tesselation *TesselStruct, class molecule *mol, int N) 683 683 { 684 if (tecplot != NULL) 685 { 686 *tecplot << "TITLE = \"3D CONVEX SHELL\"" << endl; 687 *tecplot << "VARIABLES = \"X\" \"Y\" \"Z\"" << endl; 688 *tecplot << "ZONE T=\"TRIANGLES" << N << "\", N=" 689 << TesselStruct->PointsOnBoundaryCount << ", E=" 690 << TesselStruct->TrianglesOnBoundaryCount 691 << ", DATAPACKING=POINT, ZONETYPE=FETRIANGLE" << endl; 692 int *LookupList = new int[mol->AtomCount]; 693 for (int i = 0; i < mol->AtomCount; i++) 694 LookupList[i] = -1; 695 696 // print atom coordinates 697 *out << Verbose(2) << "The following triangles were created:"; 698 int Counter = 1; 699 atom *Walker = NULL; 700 for (PointMap::iterator target = TesselStruct->PointsOnBoundary.begin(); target 701 != TesselStruct->PointsOnBoundary.end(); target++) 702 { 703 Walker = target->second->node; 704 LookupList[Walker->nr] = Counter++; 705 *tecplot << Walker->x.x[0] << " " << Walker->x.x[1] << " " 706 << Walker->x.x[2] << " " << endl; 707 } 708 *tecplot << endl; 709 // print connectivity 710 for (TriangleMap::iterator runner = 711 TesselStruct->TrianglesOnBoundary.begin(); runner 712 != TesselStruct->TrianglesOnBoundary.end(); runner++) 713 { 714 *out << " " << runner->second->endpoints[0]->node->Name << "<->" 715 << runner->second->endpoints[1]->node->Name << "<->" 716 << runner->second->endpoints[2]->node->Name; 717 *tecplot << LookupList[runner->second->endpoints[0]->node->nr] << " " 718 << LookupList[runner->second->endpoints[1]->node->nr] << " " 719 << LookupList[runner->second->endpoints[2]->node->nr] << endl; 720 } 721 delete[] (LookupList); 722 *out << endl; 723 } 684 if (tecplot != NULL) { 685 *tecplot << "TITLE = \"3D CONVEX SHELL\"" << endl; 686 *tecplot << "VARIABLES = \"X\" \"Y\" \"Z\"" << endl; 687 *tecplot << "ZONE T=\"TRIANGLES" << N << "\", N=" << TesselStruct->PointsOnBoundaryCount << ", E=" << TesselStruct->TrianglesOnBoundaryCount << ", DATAPACKING=POINT, ZONETYPE=FETRIANGLE" << endl; 688 int *LookupList = new int[mol->AtomCount]; 689 for (int i = 0; i < mol->AtomCount; i++) 690 LookupList[i] = -1; 691 692 // print atom coordinates 693 *out << Verbose(2) << "The following triangles were created:"; 694 int Counter = 1; 695 atom *Walker = NULL; 696 for (PointMap::iterator target = TesselStruct->PointsOnBoundary.begin(); target != TesselStruct->PointsOnBoundary.end(); target++) { 697 Walker = target->second->node; 698 LookupList[Walker->nr] = Counter++; 699 *tecplot << Walker->x.x[0] << " " << Walker->x.x[1] << " " << Walker->x.x[2] << " " << endl; 700 } 701 *tecplot << endl; 702 // print connectivity 703 for (TriangleMap::iterator runner = TesselStruct->TrianglesOnBoundary.begin(); runner != TesselStruct->TrianglesOnBoundary.end(); runner++) { 704 *out << " " << runner->second->endpoints[0]->node->Name << "<->" << runner->second->endpoints[1]->node->Name << "<->" << runner->second->endpoints[2]->node->Name; 705 *tecplot << LookupList[runner->second->endpoints[0]->node->nr] << " " << LookupList[runner->second->endpoints[1]->node->nr] << " " << LookupList[runner->second->endpoints[2]->node->nr] << endl; 706 } 707 delete[] (LookupList); 708 *out << endl; 709 } 724 710 } 725 711 … … 791 777 // 4. Store triangles in tecplot file 792 778 if (filename != NULL) { 793 string OutputName(filename); 794 OutputName.append(TecplotSuffix); 795 ofstream *tecplot = new ofstream(OutputName.c_str()); 796 write_tecplot_file(out, tecplot, TesselStruct, mol, 0); 797 tecplot->close(); 798 delete(tecplot); 799 ofstream *rasterplot = new ofstream(OutputName.c_str()); 800 write_raster3d_file(out, rasterplot, TesselStruct, mol); 801 rasterplot->close(); 802 delete(rasterplot); 779 if (DoTecplotOutput) { 780 string OutputName(filename); 781 OutputName.append(TecplotSuffix); 782 ofstream *tecplot = new ofstream(OutputName.c_str()); 783 write_tecplot_file(out, tecplot, TesselStruct, mol, 0); 784 tecplot->close(); 785 delete(tecplot); 786 } 787 if (DoRaster3DOutput) { 788 string OutputName(filename); 789 OutputName.append(Raster3DSuffix); 790 ofstream *rasterplot = new ofstream(OutputName.c_str()); 791 write_raster3d_file(out, rasterplot, TesselStruct, mol); 792 rasterplot->close(); 793 delete(rasterplot); 794 } 803 795 } 804 796 … … 1100 1092 << A->second->node->Name << "," 1101 1093 << baseline->second.first->second->node->Name << "," 1102 << baseline->second.second->second->node->Name << " leave "1094 << baseline->second.second->second->node->Name << " leaves " 1103 1095 << checker->second->node->Name << " outside the convex hull." 1104 1096 << endl; … … 1185 1177 * -# if the lines contains to only one triangle 1186 1178 * -# We search all points in the boundary 1187 * -# if the triangle with the baseline and the current point has the smallest of angles (comparison between normal vectors1188 1179 * -# if the triangle is in forward direction of the baseline (at most 90 degrees angle between vector orthogonal to 1189 1180 * baseline in triangle plane pointing out of the triangle and normal vector of new triangle) 1181 * -# if the triangle with the baseline and the current point has the smallest of angles (comparison between normal vectors) 1190 1182 * -# then we have a new triangle, whose baselines we again add (or increase their TriangleCount) 1191 1183 * \param *out output stream for debugging … … 1207 1199 for (LineMap::iterator baseline = LinesOnBoundary.begin(); baseline != LinesOnBoundary.end(); baseline++) 1208 1200 if (baseline->second->TrianglesCount == 1) { 1209 *out << Verbose(2) << "Current baseline is between " << *(baseline->second) << "." << endl;1210 1201 // 5a. go through each boundary point if not _both_ edges between either endpoint of the current line and this point exist (and belong to 2 triangles) 1211 1202 SmallestAngle = M_PI; 1203 1204 // get peak point with respect to this base line's only triangle 1212 1205 BTS = baseline->second->triangles.begin()->second; // there is only one triangle so far 1213 // get peak point with respect to this base line's only triangle1206 *out << Verbose(2) << "Current baseline is between " << *(baseline->second) << "." << endl; 1214 1207 for (int i = 0; i < 3; i++) 1215 1208 if ((BTS->endpoints[i] != baseline->second->endpoints[0]) && (BTS->endpoints[i] != baseline->second->endpoints[1])) 1216 1209 peak = BTS->endpoints[i]; 1217 1210 *out << Verbose(3) << " and has peak " << *peak << "." << endl; 1211 1212 // prepare some auxiliary vectors 1213 Vector BaseLineCenter, BaseLine; 1214 BaseLineCenter.CopyVector(&baseline->second->endpoints[0]->node->x); 1215 BaseLineCenter.AddVector(&baseline->second->endpoints[0]->node->x); 1216 BaseLineCenter.Scale(1. / 2.); // points now to center of base line 1217 BaseLine.CopyVector(&baseline->second->endpoints[0]->node->x); 1218 BaseLine.SubtractVector(&baseline->second->endpoints[1]->node->x); 1219 1218 1220 // offset to center of triangle 1219 1221 CenterVector.Zero(); … … 1222 1224 CenterVector.Scale(1. / 3.); 1223 1225 *out << Verbose(4) << "CenterVector of base triangle is " << CenterVector << endl; 1226 1227 // normal vector of triangle 1224 1228 NormalVector.CopyVector(MolCenter); 1225 1229 NormalVector.SubtractVector(&CenterVector); 1226 // normal vector of triangle1227 1230 BTS->GetNormalVector(NormalVector); 1231 NormalVector.CopyVector(&BTS->NormalVector); 1228 1232 *out << Verbose(4) << "NormalVector of base triangle is " << NormalVector << endl; 1233 1229 1234 // vector in propagation direction (out of triangle) 1230 1235 // project center vector onto triangle plane (points from intersection plane-NormalVector to plane-CenterVector intersection) 1231 TempVector.CopyVector(&baseline->second->endpoints[0]->node->x); 1232 TempVector.SubtractVector(&baseline->second->endpoints[1]->node->x); 1233 PropagationVector.MakeNormalVector(&TempVector, &NormalVector); 1236 PropagationVector.MakeNormalVector(&BaseLine, &NormalVector); 1234 1237 TempVector.CopyVector(&CenterVector); 1235 1238 TempVector.SubtractVector(&baseline->second->endpoints[0]->node->x); // TempVector is vector on triangle plane pointing from one baseline egde towards center! … … 1237 1240 if (PropagationVector.Projection(&TempVector) > 0) // make sure normal propagation vector points outward from baseline 1238 1241 PropagationVector.Scale(-1.); 1239 *out << Verbose(4) << "PropagationVector of base triangle is "; 1240 PropagationVector.Output(out); 1241 *out << endl; 1242 *out << Verbose(4) << "PropagationVector of base triangle is " << PropagationVector << endl; 1242 1243 winner = PointsOnBoundary.end(); 1243 for (PointMap::iterator target = PointsOnBoundary.begin(); target 1244 != PointsOnBoundary.end(); target++) 1245 if ((target->second != baseline->second->endpoints[0]) 1246 && (target->second != baseline->second->endpoints[1])) 1247 { // don't take the same endpoints 1248 *out << Verbose(3) << "Target point is " << *(target->second) 1249 << ":"; 1250 bool continueflag = true; 1251 1252 VirtualNormalVector.CopyVector( 1253 &baseline->second->endpoints[0]->node->x); 1254 VirtualNormalVector.AddVector( 1255 &baseline->second->endpoints[0]->node->x); 1256 VirtualNormalVector.Scale(-1. / 2.); // points now to center of base line 1257 VirtualNormalVector.AddVector(&target->second->node->x); // points from center of base line to target 1258 TempAngle = VirtualNormalVector.Angle(&PropagationVector); 1259 continueflag = continueflag && (TempAngle < (M_PI/2.)); // no bends bigger than Pi/2 (90 degrees) 1260 if (!continueflag) 1261 { 1262 *out << Verbose(4) 1263 << "Angle between propagation direction and base line to " 1264 << *(target->second) << " is " << TempAngle 1265 << ", bad direction!" << endl; 1266 continue; 1267 } 1268 else 1269 *out << Verbose(4) 1270 << "Angle between propagation direction and base line to " 1271 << *(target->second) << " is " << TempAngle 1272 << ", good direction!" << endl; 1273 LineChecker[0] = baseline->second->endpoints[0]->lines.find( 1274 target->first); 1275 LineChecker[1] = baseline->second->endpoints[1]->lines.find( 1276 target->first); 1277 // if (LineChecker[0] != baseline->second->endpoints[0]->lines.end()) 1278 // *out << Verbose(4) << *(baseline->second->endpoints[0]) << " has line " << *(LineChecker[0]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[0]->second->TrianglesCount << " triangles." << endl; 1279 // else 1280 // *out << Verbose(4) << *(baseline->second->endpoints[0]) << " has no line to " << *(target->second) << " as endpoint." << endl; 1281 // if (LineChecker[1] != baseline->second->endpoints[1]->lines.end()) 1282 // *out << Verbose(4) << *(baseline->second->endpoints[1]) << " has line " << *(LineChecker[1]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[1]->second->TrianglesCount << " triangles." << endl; 1283 // else 1284 // *out << Verbose(4) << *(baseline->second->endpoints[1]) << " has no line to " << *(target->second) << " as endpoint." << endl; 1285 // check first endpoint (if any connecting line goes to target or at least not more than 1) 1286 continueflag = continueflag && (((LineChecker[0] 1287 == baseline->second->endpoints[0]->lines.end()) 1288 || (LineChecker[0]->second->TrianglesCount == 1))); 1289 if (!continueflag) 1290 { 1291 *out << Verbose(4) << *(baseline->second->endpoints[0]) 1292 << " has line " << *(LineChecker[0]->second) 1293 << " to " << *(target->second) 1294 << " as endpoint with " 1295 << LineChecker[0]->second->TrianglesCount 1296 << " triangles." << endl; 1297 continue; 1298 } 1299 // check second endpoint (if any connecting line goes to target or at least not more than 1) 1300 continueflag = continueflag && (((LineChecker[1] 1301 == baseline->second->endpoints[1]->lines.end()) 1302 || (LineChecker[1]->second->TrianglesCount == 1))); 1303 if (!continueflag) 1304 { 1305 *out << Verbose(4) << *(baseline->second->endpoints[1]) 1306 << " has line " << *(LineChecker[1]->second) 1307 << " to " << *(target->second) 1308 << " as endpoint with " 1309 << LineChecker[1]->second->TrianglesCount 1310 << " triangles." << endl; 1311 continue; 1312 } 1313 // check whether the envisaged triangle does not already exist (if both lines exist and have same endpoint) 1314 continueflag = continueflag && (!(((LineChecker[0] 1315 != baseline->second->endpoints[0]->lines.end()) 1316 && (LineChecker[1] 1317 != baseline->second->endpoints[1]->lines.end()) 1318 && (GetCommonEndpoint(LineChecker[0]->second, 1319 LineChecker[1]->second) == peak)))); 1320 if (!continueflag) 1321 { 1322 *out << Verbose(4) << "Current target is peak!" << endl; 1323 continue; 1324 } 1325 // in case NOT both were found 1326 if (continueflag) 1327 { // create virtually this triangle, get its normal vector, calculate angle 1328 flag = true; 1329 VirtualNormalVector.MakeNormalVector( 1330 &baseline->second->endpoints[0]->node->x, 1331 &baseline->second->endpoints[1]->node->x, 1332 &target->second->node->x); 1333 // make it always point inward 1334 if (baseline->second->endpoints[0]->node->x.Projection( 1335 &VirtualNormalVector) > 0) 1336 VirtualNormalVector.Scale(-1.); 1337 // calculate angle 1338 TempAngle = NormalVector.Angle(&VirtualNormalVector); 1339 *out << Verbose(4) << "NormalVector is "; 1340 VirtualNormalVector.Output(out); 1341 *out << " and the angle is " << TempAngle << "." << endl; 1342 if (SmallestAngle > TempAngle) 1343 { // set to new possible winner 1344 SmallestAngle = TempAngle; 1345 winner = target; 1346 } 1347 } 1244 1245 // loop over all points and calculate angle between normal vector of new and present triangle 1246 for (PointMap::iterator target = PointsOnBoundary.begin(); target != PointsOnBoundary.end(); target++) { 1247 if ((target->second != baseline->second->endpoints[0]) && (target->second != baseline->second->endpoints[1])) { // don't take the same endpoints 1248 *out << Verbose(3) << "Target point is " << *(target->second) << ":" << endl; 1249 1250 // first check direction, so that triangles don't intersect 1251 VirtualNormalVector.CopyVector(&target->second->node->x); 1252 VirtualNormalVector.SubtractVector(&BaseLineCenter); // points from center of base line to target 1253 VirtualNormalVector.ProjectOntoPlane(&NormalVector); 1254 TempAngle = VirtualNormalVector.Angle(&PropagationVector); 1255 if (TempAngle > (M_PI/2.)) { // no bends bigger than Pi/2 (90 degrees) 1256 *out << Verbose(4) << "Angle on triangle plane between propagation direction and base line to " << *(target->second) << " is " << TempAngle << ", bad direction!" << endl; 1257 continue; 1258 } else 1259 *out << Verbose(4) << "Angle on triangle plane between propagation direction and base line to " << *(target->second) << " is " << TempAngle << ", good direction!" << endl; 1260 1261 // check first and second endpoint (if any connecting line goes to target has at least not more than 1 triangle) 1262 LineChecker[0] = baseline->second->endpoints[0]->lines.find(target->first); 1263 LineChecker[1] = baseline->second->endpoints[1]->lines.find(target->first); 1264 if (((LineChecker[0] != baseline->second->endpoints[0]->lines.end()) && (LineChecker[0]->second->TrianglesCount == 2))) { 1265 *out << Verbose(4) << *(baseline->second->endpoints[0]) << " has line " << *(LineChecker[0]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[0]->second->TrianglesCount << " triangles." << endl; 1266 continue; 1348 1267 } 1268 if (((LineChecker[1] != baseline->second->endpoints[1]->lines.end()) && (LineChecker[1]->second->TrianglesCount == 2))) { 1269 *out << Verbose(4) << *(baseline->second->endpoints[1]) << " has line " << *(LineChecker[1]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[1]->second->TrianglesCount << " triangles." << endl; 1270 continue; 1271 } 1272 1273 // check whether the envisaged triangle does not already exist (if both lines exist and have same endpoint) 1274 if ((((LineChecker[0] != baseline->second->endpoints[0]->lines.end()) && (LineChecker[1] != baseline->second->endpoints[1]->lines.end()) && (GetCommonEndpoint(LineChecker[0]->second, LineChecker[1]->second) == peak)))) { 1275 *out << Verbose(4) << "Current target is peak!" << endl; 1276 continue; 1277 } 1278 1279 // in case NOT both were found, create virtually this triangle, get its normal vector, calculate angle 1280 flag = true; 1281 VirtualNormalVector.MakeNormalVector(&baseline->second->endpoints[0]->node->x, &baseline->second->endpoints[1]->node->x, &target->second->node->x); 1282 TempVector.CopyVector(&baseline->second->endpoints[0]->node->x); 1283 TempVector.AddVector(&baseline->second->endpoints[1]->node->x); 1284 TempVector.AddVector(&target->second->node->x); 1285 TempVector.Scale(1./3.); 1286 TempVector.SubtractVector(MolCenter); 1287 // make it always point outward 1288 if (VirtualNormalVector.Projection(&TempVector) < 0) 1289 VirtualNormalVector.Scale(-1.); 1290 // calculate angle 1291 TempAngle = NormalVector.Angle(&VirtualNormalVector); 1292 *out << Verbose(4) << "NormalVector is " << VirtualNormalVector << " and the angle is " << TempAngle << "." << endl; 1293 if (SmallestAngle > TempAngle) { // set to new possible winner 1294 SmallestAngle = TempAngle; 1295 winner = target; 1296 } 1297 } 1298 } // end of loop over all boundary points 1299 1349 1300 // 5b. The point of the above whose triangle has the greatest angle with the triangle the current line belongs to (it only belongs to one, remember!): New triangle 1350 if (winner != PointsOnBoundary.end()) 1351 { 1352 *out << Verbose(2) << "Winning target point is " 1353 << *(winner->second) << " with angle " << SmallestAngle 1354 << "." << endl; 1355 // create the lins of not yet present 1356 BLS[0] = baseline->second; 1357 // 5c. add lines to the line set if those were new (not yet part of a triangle), delete lines that belong to two triangles) 1358 LineChecker[0] = baseline->second->endpoints[0]->lines.find( 1359 winner->first); 1360 LineChecker[1] = baseline->second->endpoints[1]->lines.find( 1361 winner->first); 1362 if (LineChecker[0] 1363 == baseline->second->endpoints[0]->lines.end()) 1364 { // create 1365 BPS[0] = baseline->second->endpoints[0]; 1366 BPS[1] = winner->second; 1367 BLS[1] = new class BoundaryLineSet(BPS, 1368 LinesOnBoundaryCount); 1369 LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, 1370 BLS[1])); 1371 LinesOnBoundaryCount++; 1372 } 1373 else 1374 BLS[1] = LineChecker[0]->second; 1375 if (LineChecker[1] 1376 == baseline->second->endpoints[1]->lines.end()) 1377 { // create 1378 BPS[0] = baseline->second->endpoints[1]; 1379 BPS[1] = winner->second; 1380 BLS[2] = new class BoundaryLineSet(BPS, 1381 LinesOnBoundaryCount); 1382 LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, 1383 BLS[2])); 1384 LinesOnBoundaryCount++; 1385 } 1386 else 1387 BLS[2] = LineChecker[1]->second; 1388 BTS = new class BoundaryTriangleSet(BLS, 1389 TrianglesOnBoundaryCount); 1390 TrianglesOnBoundary.insert(TrianglePair( 1391 TrianglesOnBoundaryCount, BTS)); 1392 TrianglesOnBoundaryCount++; 1393 } 1394 else 1395 { 1396 *out << Verbose(1) 1397 << "I could not determine a winner for this baseline " 1398 << *(baseline->second) << "." << endl; 1399 } 1301 if (winner != PointsOnBoundary.end()) { 1302 *out << Verbose(2) << "Winning target point is " << *(winner->second) << " with angle " << SmallestAngle << "." << endl; 1303 // create the lins of not yet present 1304 BLS[0] = baseline->second; 1305 // 5c. add lines to the line set if those were new (not yet part of a triangle), delete lines that belong to two triangles) 1306 LineChecker[0] = baseline->second->endpoints[0]->lines.find(winner->first); 1307 LineChecker[1] = baseline->second->endpoints[1]->lines.find(winner->first); 1308 if (LineChecker[0] == baseline->second->endpoints[0]->lines.end()) { // create 1309 BPS[0] = baseline->second->endpoints[0]; 1310 BPS[1] = winner->second; 1311 BLS[1] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount); 1312 LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BLS[1])); 1313 LinesOnBoundaryCount++; 1314 } else 1315 BLS[1] = LineChecker[0]->second; 1316 if (LineChecker[1] == baseline->second->endpoints[1]->lines.end()) { // create 1317 BPS[0] = baseline->second->endpoints[1]; 1318 BPS[1] = winner->second; 1319 BLS[2] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount); 1320 LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BLS[2])); 1321 LinesOnBoundaryCount++; 1322 } else 1323 BLS[2] = LineChecker[1]->second; 1324 BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount); 1325 TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS)); 1326 TrianglesOnBoundaryCount++; 1327 } else { 1328 *out << Verbose(1) << "I could not determine a winner for this baseline " << *(baseline->second) << "." << endl; 1329 } 1400 1330 1401 1331 // 5d. If the set of lines is not yet empty, go to 5. and continue 1402 } 1403 else 1404 *out << Verbose(2) << "Baseline candidate " << *(baseline->second) 1405 << " has a triangle count of " 1406 << baseline->second->TrianglesCount << "." << endl; 1332 } else 1333 *out << Verbose(2) << "Baseline candidate " << *(baseline->second) << " has a triangle count of " << baseline->second->TrianglesCount << "." << endl; 1407 1334 } while (flag); 1408 1335 delete(MolCenter); … … 2795 2722 } 2796 2723 } 2797 if ( 1) { //failflag) {2724 if (filename != 0) { 2798 2725 *out << Verbose(1) << "Writing final tecplot file\n"; 2799 2726 if (DoTecplotOutput) { -
src/builder.cpp
rfcc9b0 ra567c3 1813 1813 LinkedCell LCList(mol, 10.); 1814 1814 class Tesselation *TesselStruct = NULL; 1815 Find_convex_border((ofstream *)&cout, mol, TesselStruct, &LCList, NULL);1815 Find_convex_border((ofstream *)&cout, mol, TesselStruct, &LCList, argv[argptr]); 1816 1816 double clustervolume = VolumeOfConvexEnvelope((ofstream *)&cout, TesselStruct, &configuration); 1817 1817 delete(TesselStruct);
Note:
See TracChangeset
for help on using the changeset viewer.