Changeset 042f82 for src/boundary.cpp
- Timestamp:
- Jul 23, 2009, 2:21:57 PM (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:
- 36ec71
- Parents:
- 205ccd
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
src/boundary.cpp
r205ccd r042f82 239 239 240 240 CandidateForTesselation::CandidateForTesselation( 241 241 atom *candidate, BoundaryLineSet* line, Vector OptCandidateCenter, Vector OtherOptCandidateCenter 242 242 ) { 243 244 245 246 243 point = candidate; 244 BaseLine = line; 245 OptCenter.CopyVector(&OptCandidateCenter); 246 OtherOptCenter.CopyVector(&OtherOptCandidateCenter); 247 247 } 248 248 249 249 CandidateForTesselation::~CandidateForTesselation() { 250 251 250 point = NULL; 251 BaseLine = NULL; 252 252 } 253 253 … … 1904 1904 Vector CirclePlaneNormal; // normal vector defining the plane this circle lives in 1905 1905 Vector SphereCenter; 1906 Vector NewSphereCenter; 1907 Vector OtherNewSphereCenter; 1906 Vector NewSphereCenter; // center of the sphere defined by the two points of BaseLine and the one of Candidate, first possibility 1907 Vector OtherNewSphereCenter; // center of the sphere defined by the two points of BaseLine and the one of Candidate, second possibility 1908 1908 Vector NewNormalVector; // normal vector of the Candidate's triangle 1909 1909 Vector helper, OptCandidateCenter, OtherOptCandidateCenter; … … 1985 1985 1986 1986 if ((NewNormalVector.MakeNormalVector(&(BaseLine->endpoints[0]->node->x), &(BaseLine->endpoints[1]->node->x), &(Candidate->x))) 1987 1987 && (fabs(NewNormalVector.ScalarProduct(&NewNormalVector)) > HULLEPSILON) 1988 1988 ) { 1989 1989 helper.CopyVector(&NewNormalVector); … … 2072 2072 //cout << Verbose(2) << "INFO: Sorting candidate list ..." << endl; 2073 2073 if (candidates->size() > 1) { 2074 2075 2074 candidates->unique(); 2075 candidates->sort(sortCandidates); 2076 2076 } 2077 2077 … … 2080 2080 2081 2081 struct Intersection { 2082 2083 2084 2085 2082 Vector x1; 2083 Vector x2; 2084 Vector x3; 2085 Vector x4; 2086 2086 }; 2087 2087 … … 2093 2093 */ 2094 2094 double MinIntersectDistance(const gsl_vector * x, void *params) { 2095 2096 2097 2098 2099 2100 2101 2102 2103 2104 2105 2106 2107 2108 2109 2110 2111 2112 2113 2114 2115 2116 2117 2095 double retval = 0; 2096 struct Intersection *I = (struct Intersection *)params; 2097 Vector intersection; 2098 Vector SideA,SideB,HeightA, HeightB; 2099 for (int i=0;i<NDIM;i++) 2100 intersection.x[i] = gsl_vector_get(x, i); 2101 2102 SideA.CopyVector(&(I->x1)); 2103 SideA.SubtractVector(&I->x2); 2104 HeightA.CopyVector(&intersection); 2105 HeightA.SubtractVector(&I->x1); 2106 HeightA.ProjectOntoPlane(&SideA); 2107 2108 SideB.CopyVector(&I->x3); 2109 SideB.SubtractVector(&I->x4); 2110 HeightB.CopyVector(&intersection); 2111 HeightB.SubtractVector(&I->x3); 2112 HeightB.ProjectOntoPlane(&SideB); 2113 2114 retval = HeightA.ScalarProduct(&HeightA) + HeightB.ScalarProduct(&HeightB); 2115 //cout << Verbose(2) << "MinIntersectDistance called, result: " << retval << endl; 2116 2117 return retval; 2118 2118 }; 2119 2119 … … 2132 2132 */ 2133 2133 bool existsIntersection(Vector point1, Vector point2, Vector point3, Vector point4) { 2134 2135 2136 2137 2138 2139 2140 2134 bool result; 2135 2136 struct Intersection par; 2137 par.x1.CopyVector(&point1); 2138 par.x2.CopyVector(&point2); 2139 par.x3.CopyVector(&point3); 2140 par.x4.CopyVector(&point4); 2141 2141 2142 2142 const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex; … … 2179 2179 2180 2180 if (status == GSL_SUCCESS) { 2181 2181 cout << Verbose(2) << "converged to minimum" << endl; 2182 2182 } 2183 2183 } while (status == GSL_CONTINUE && iter < 100); 2184 2184 2185 2185 // check whether intersection is in between or not 2186 2187 2188 2189 2190 2191 2192 2193 2194 2195 2196 2197 2198 2199 2200 2201 2202 2203 2204 2205 2206 2207 2208 2209 2210 2211 2212 2213 2214 2215 2216 2217 2218 2186 Vector intersection, SideA, SideB, HeightA, HeightB; 2187 double t1, t2; 2188 for (int i = 0; i < NDIM; i++) { 2189 intersection.x[i] = gsl_vector_get(s->x, i); 2190 } 2191 2192 SideA.CopyVector(&par.x2); 2193 SideA.SubtractVector(&par.x1); 2194 HeightA.CopyVector(&intersection); 2195 HeightA.SubtractVector(&par.x1); 2196 2197 t1 = HeightA.Projection(&SideA)/SideA.ScalarProduct(&SideA); 2198 2199 SideB.CopyVector(&par.x4); 2200 SideB.SubtractVector(&par.x3); 2201 HeightB.CopyVector(&intersection); 2202 HeightB.SubtractVector(&par.x3); 2203 2204 t2 = HeightB.Projection(&SideB)/SideB.ScalarProduct(&SideB); 2205 2206 cout << Verbose(2) << "Intersection " << intersection << " is at " 2207 << t1 << " for (" << point1 << "," << point2 << ") and at " 2208 << t2 << " for (" << point3 << "," << point4 << "): "; 2209 2210 if (((t1 >= 0) && (t1 <= 1)) && ((t2 >= 0) && (t2 <= 1))) { 2211 cout << "true intersection." << endl; 2212 result = true; 2213 } else { 2214 cout << "intersection out of region of interest." << endl; 2215 result = false; 2216 } 2217 2218 // free minimizer stuff 2219 2219 gsl_vector_free(x); 2220 2220 gsl_vector_free(ss); 2221 2221 gsl_multimin_fminimizer_free(s); 2222 2222 2223 2223 return result; 2224 2224 } 2225 2225 … … 2603 2603 cout << Verbose(1) << "Third Points are "; 2604 2604 for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) { 2605 2605 cout << " " << *(*it)->point; 2606 2606 } 2607 2607 cout << endl; … … 2609 2609 BoundaryLineSet *BaseRay = &Line; 2610 2610 for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) { 2611 2612 2613 2614 2615 2616 2617 2618 2619 2620 2621 2622 2623 2624 2611 cout << Verbose(1) << " Third point candidate is " << *(*it)->point 2612 << " with circumsphere's center at " << (*it)->OptCenter << "." << endl; 2613 cout << Verbose(1) << " Baseline is " << *BaseRay << endl; 2614 2615 // check whether all edges of the new triangle still have space for one more triangle (i.e. TriangleCount <2) 2616 atom *AtomCandidates[3]; 2617 AtomCandidates[0] = (*it)->point; 2618 AtomCandidates[1] = BaseRay->endpoints[0]->node; 2619 AtomCandidates[2] = BaseRay->endpoints[1]->node; 2620 int existentTrianglesCount = CheckPresenceOfTriangle(out, AtomCandidates); 2621 2622 BTS = NULL; 2623 // If there is no triangle, add it regularly. 2624 if (existentTrianglesCount == 0) { 2625 2625 AddTrianglePoint((*it)->point, 0); 2626 2626 AddTrianglePoint(BaseRay->endpoints[0]->node, 1); … … 2640 2640 << " for this triangle ... " << endl; 2641 2641 //cout << Verbose(1) << "We have "<< TrianglesOnBoundaryCount << " for line " << *BaseRay << "." << endl; 2642 2642 } else if (existentTrianglesCount == 1) { // If there is a planar region within the structure, we need this triangle a second time. 2643 2643 AddTrianglePoint((*it)->point, 0); 2644 2644 AddTrianglePoint(BaseRay->endpoints[0]->node, 1); … … 2679 2679 } 2680 2680 2681 2681 if ((result) && (existentTrianglesCount < 2) && (DoSingleStepOutput && (TrianglesOnBoundaryCount % 1 == 0))) { // if we have a new triangle and want to output each new triangle configuration 2682 2682 sprintf(NumberName, "-%04d-%s_%s_%s", TriangleFilesWritten, BTS->endpoints[0]->node->Name, BTS->endpoints[1]->node->Name, BTS->endpoints[2]->node->Name); 2683 2683 if (DoTecplotOutput) { … … 2726 2726 if (DoTecplotOutput || DoRaster3DOutput) 2727 2727 TriangleFilesWritten++; 2728 2728 } 2729 2729 2730 2730 // set baseline to new ray from ref point (here endpoints[0]->node) to current candidate (here (*it)->point)) 2731 2731 BaseRay = BLS[0]; 2732 2732 } 2733 2733 … … 2747 2747 */ 2748 2748 bool sortCandidates(CandidateForTesselation* candidate1, CandidateForTesselation* candidate2) { 2749 2750 2751 2752 2753 2754 2755 2756 2757 2758 2749 Vector BaseLineVector, OrthogonalVector, helper; 2750 if (candidate1->BaseLine != candidate2->BaseLine) { // sanity check 2751 cout << Verbose(0) << "ERROR: sortCandidates was called for two different baselines: " << candidate1->BaseLine << " and " << candidate2->BaseLine << "." << endl; 2752 //return false; 2753 exit(1); 2754 } 2755 // create baseline vector 2756 BaseLineVector.CopyVector(&(candidate1->BaseLine->endpoints[1]->node->x)); 2757 BaseLineVector.SubtractVector(&(candidate1->BaseLine->endpoints[0]->node->x)); 2758 BaseLineVector.Normalize(); 2759 2759 2760 2760 // create normal in-plane vector to cope with acos() non-uniqueness on [0,2pi] (note that is pointing in the "right" direction already, hence ">0" test!) 2761 2762 2763 2764 2765 2766 2761 helper.CopyVector(&(candidate1->BaseLine->endpoints[0]->node->x)); 2762 helper.SubtractVector(&(candidate1->point->x)); 2763 OrthogonalVector.CopyVector(&helper); 2764 helper.VectorProduct(&BaseLineVector); 2765 OrthogonalVector.SubtractVector(&helper); 2766 OrthogonalVector.Normalize(); 2767 2767 2768 2768 // calculate both angles and correct with in-plane vector 2769 2770 2771 2769 helper.CopyVector(&(candidate1->point->x)); 2770 helper.SubtractVector(&(candidate1->BaseLine->endpoints[0]->node->x)); 2771 double phi = BaseLineVector.Angle(&helper); 2772 2772 if (OrthogonalVector.ScalarProduct(&helper) > 0) { 2773 2773 phi = 2.*M_PI - phi; 2774 2774 } 2775 2775 helper.CopyVector(&(candidate2->point->x)); 2776 2777 2776 helper.SubtractVector(&(candidate1->BaseLine->endpoints[0]->node->x)); 2777 double psi = BaseLineVector.Angle(&helper); 2778 2778 if (OrthogonalVector.ScalarProduct(&helper) > 0) { 2779 2780 2781 2782 2783 2784 2785 2786 2779 psi = 2.*M_PI - psi; 2780 } 2781 2782 cout << Verbose(2) << *candidate1->point << " has angle " << phi << endl; 2783 cout << Verbose(2) << *candidate2->point << " has angle " << psi << endl; 2784 2785 // return comparison 2786 return phi < psi; 2787 2787 } 2788 2788
Note:
See TracChangeset
for help on using the changeset viewer.