Changeset 042f82 for src/boundary.cpp


Ignore:
Timestamp:
Jul 23, 2009, 2:21:57 PM (16 years ago)
Author:
Frederik Heber <heber@…>
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
Message:

fixed indentation from tabs to two spaces.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • src/boundary.cpp

    r205ccd r042f82  
    239239
    240240CandidateForTesselation::CandidateForTesselation(
    241         atom *candidate, BoundaryLineSet* line, Vector OptCandidateCenter, Vector OtherOptCandidateCenter
     241  atom *candidate, BoundaryLineSet* line, Vector OptCandidateCenter, Vector OtherOptCandidateCenter
    242242) {
    243         point = candidate;
    244         BaseLine = line;
    245         OptCenter.CopyVector(&OptCandidateCenter);
    246         OtherOptCenter.CopyVector(&OtherOptCandidateCenter);
     243  point = candidate;
     244  BaseLine = line;
     245  OptCenter.CopyVector(&OptCandidateCenter);
     246  OtherOptCenter.CopyVector(&OtherOptCandidateCenter);
    247247}
    248248
    249249CandidateForTesselation::~CandidateForTesselation() {
    250         point = NULL;
    251         BaseLine = NULL;
     250  point = NULL;
     251  BaseLine = NULL;
    252252}
    253253
     
    19041904  Vector CirclePlaneNormal; // normal vector defining the plane this circle lives in
    19051905  Vector SphereCenter;
    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
     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
    19081908  Vector NewNormalVector;   // normal vector of the Candidate's triangle
    19091909  Vector helper, OptCandidateCenter, OtherOptCandidateCenter;
     
    19851985
    19861986                  if ((NewNormalVector.MakeNormalVector(&(BaseLine->endpoints[0]->node->x), &(BaseLine->endpoints[1]->node->x), &(Candidate->x)))
    1987                         && (fabs(NewNormalVector.ScalarProduct(&NewNormalVector)) > HULLEPSILON)
     1987                  && (fabs(NewNormalVector.ScalarProduct(&NewNormalVector)) > HULLEPSILON)
    19881988                  ) {
    19891989                    helper.CopyVector(&NewNormalVector);
     
    20722072  //cout << Verbose(2) << "INFO: Sorting candidate list ..." << endl;
    20732073  if (candidates->size() > 1) {
    2074           candidates->unique();
    2075           candidates->sort(sortCandidates);
     2074    candidates->unique();
     2075    candidates->sort(sortCandidates);
    20762076  }
    20772077
     
    20802080
    20812081struct Intersection {
    2082         Vector x1;
    2083         Vector x2;
    2084         Vector x3;
    2085         Vector x4;
     2082  Vector x1;
     2083  Vector x2;
     2084  Vector x3;
     2085  Vector x4;
    20862086};
    20872087
     
    20932093 */
    20942094double MinIntersectDistance(const gsl_vector * x, void *params) {
    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;
     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;
    21182118};
    21192119
     
    21322132 */
    21332133bool existsIntersection(Vector point1, Vector point2, Vector point3, Vector point4) {
    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);
     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);
    21412141
    21422142    const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex;
     
    21792179
    21802180        if (status == GSL_SUCCESS) {
    2181                 cout << Verbose(2) << "converged to minimum" <<  endl;
     2181          cout << Verbose(2) << "converged to minimum" <<  endl;
    21822182        }
    21832183    } while (status == GSL_CONTINUE && iter < 100);
    21842184
    21852185    // check whether intersection is in between or not
    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
     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
    22192219    gsl_vector_free(x);
    22202220    gsl_vector_free(ss);
    22212221    gsl_multimin_fminimizer_free(s);
    22222222
    2223         return result;
     2223  return result;
    22242224}
    22252225
     
    26032603  cout << Verbose(1) << "Third Points are ";
    26042604  for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
    2605           cout << " " << *(*it)->point;
     2605    cout << " " << *(*it)->point;
    26062606  }
    26072607  cout << endl;
     
    26092609  BoundaryLineSet *BaseRay = &Line;
    26102610  for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
    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) {
     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) {
    26252625      AddTrianglePoint((*it)->point, 0);
    26262626      AddTrianglePoint(BaseRay->endpoints[0]->node, 1);
     
    26402640        << " for this triangle ... " << endl;
    26412641      //cout << Verbose(1) << "We have "<< TrianglesOnBoundaryCount << " for line " << *BaseRay << "." << endl;
    2642           } else if (existentTrianglesCount == 1) { // If there is a planar region within the structure, we need this triangle a second time.
     2642    } else if (existentTrianglesCount == 1) { // If there is a planar region within the structure, we need this triangle a second time.
    26432643        AddTrianglePoint((*it)->point, 0);
    26442644        AddTrianglePoint(BaseRay->endpoints[0]->node, 1);
     
    26792679    }
    26802680
    2681           if ((result) && (existentTrianglesCount < 2) && (DoSingleStepOutput && (TrianglesOnBoundaryCount % 1 == 0))) { // if we have a new triangle and want to output each new triangle configuration
     2681    if ((result) && (existentTrianglesCount < 2) && (DoSingleStepOutput && (TrianglesOnBoundaryCount % 1 == 0))) { // if we have a new triangle and want to output each new triangle configuration
    26822682      sprintf(NumberName, "-%04d-%s_%s_%s", TriangleFilesWritten, BTS->endpoints[0]->node->Name, BTS->endpoints[1]->node->Name, BTS->endpoints[2]->node->Name);
    26832683      if (DoTecplotOutput) {
     
    27262726      if (DoTecplotOutput || DoRaster3DOutput)
    27272727        TriangleFilesWritten++;
    2728           }
     2728    }
    27292729
    27302730    // set baseline to new ray from ref point (here endpoints[0]->node) to current candidate (here (*it)->point))
    2731           BaseRay = BLS[0];
     2731    BaseRay = BLS[0];
    27322732  }
    27332733
     
    27472747 */
    27482748bool sortCandidates(CandidateForTesselation* candidate1, CandidateForTesselation* candidate2) {
    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();
     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();
    27592759
    27602760  // 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         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();
     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();
    27672767
    27682768  // calculate both angles and correct with in-plane vector
    2769         helper.CopyVector(&(candidate1->point->x));
    2770         helper.SubtractVector(&(candidate1->BaseLine->endpoints[0]->node->x));
    2771         double phi = BaseLineVector.Angle(&helper);
     2769  helper.CopyVector(&(candidate1->point->x));
     2770  helper.SubtractVector(&(candidate1->BaseLine->endpoints[0]->node->x));
     2771  double phi = BaseLineVector.Angle(&helper);
    27722772  if (OrthogonalVector.ScalarProduct(&helper) > 0) {
    27732773    phi = 2.*M_PI - phi;
    27742774  }
    27752775  helper.CopyVector(&(candidate2->point->x));
    2776         helper.SubtractVector(&(candidate1->BaseLine->endpoints[0]->node->x));
    2777         double psi = BaseLineVector.Angle(&helper);
     2776  helper.SubtractVector(&(candidate1->BaseLine->endpoints[0]->node->x));
     2777  double psi = BaseLineVector.Angle(&helper);
    27782778  if (OrthogonalVector.ScalarProduct(&helper) > 0) {
    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;
     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;
    27872787}
    27882788
Note: See TracChangeset for help on using the changeset viewer.