Ignore:
Timestamp:
Apr 29, 2010, 1:55:21 PM (15 years ago)
Author:
Tillmann Crueger <crueger@…>
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:
d79639
Parents:
632bc3 (diff), 753f02 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge branch 'VectorRefactoring' into StructureRefactoring

Conflicts:

molecuilder/src/Legacy/oldmenu.cpp
molecuilder/src/Makefile.am
molecuilder/src/analysis_correlation.cpp
molecuilder/src/boundary.cpp
molecuilder/src/builder.cpp
molecuilder/src/config.cpp
molecuilder/src/ellipsoid.cpp
molecuilder/src/linkedcell.cpp
molecuilder/src/molecule.cpp
molecuilder/src/molecule_fragmentation.cpp
molecuilder/src/molecule_geometry.cpp
molecuilder/src/molecule_graph.cpp
molecuilder/src/moleculelist.cpp
molecuilder/src/tesselation.cpp
molecuilder/src/tesselationhelpers.cpp
molecuilder/src/unittests/AnalysisCorrelationToSurfaceUnitTest.cpp
molecuilder/src/unittests/bondgraphunittest.cpp
molecuilder/src/vector.cpp
molecuilder/src/vector.hpp

File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/tesselationhelpers.cpp

    r632bc3 r8cbb97  
    1414#include "tesselationhelpers.hpp"
    1515#include "vector.hpp"
     16#include "vector_ops.hpp"
    1617#include "verbose.hpp"
    1718
    … …  
    5354
    5455  for(int i=0;i<3;i++) {
    55     gsl_matrix_set(A, i, 0, a.x[i]);
    56     gsl_matrix_set(A, i, 1, b.x[i]);
    57     gsl_matrix_set(A, i, 2, c.x[i]);
     56    gsl_matrix_set(A, i, 0, a[i]);
     57    gsl_matrix_set(A, i, 1, b[i]);
     58    gsl_matrix_set(A, i, 2, c[i]);
    5859  }
    5960  m11 = DetGet(A, 1);
    6061
    6162  for(int i=0;i<3;i++) {
    62     gsl_matrix_set(A, i, 0, a.x[i]*a.x[i] + b.x[i]*b.x[i] + c.x[i]*c.x[i]);
    63     gsl_matrix_set(A, i, 1, b.x[i]);
    64     gsl_matrix_set(A, i, 2, c.x[i]);
     63    gsl_matrix_set(A, i, 0, a[i]*a[i] + b[i]*b[i] + c[i]*c[i]);
     64    gsl_matrix_set(A, i, 1, b[i]);
     65    gsl_matrix_set(A, i, 2, c[i]);
    6566  }
    6667  m12 = DetGet(A, 1);
    6768
    6869  for(int i=0;i<3;i++) {
    69     gsl_matrix_set(A, i, 0, a.x[i]*a.x[i] + b.x[i]*b.x[i] + c.x[i]*c.x[i]);
    70     gsl_matrix_set(A, i, 1, a.x[i]);
    71     gsl_matrix_set(A, i, 2, c.x[i]);
     70    gsl_matrix_set(A, i, 0, a[i]*a[i] + b[i]*b[i] + c[i]*c[i]);
     71    gsl_matrix_set(A, i, 1, a[i]);
     72    gsl_matrix_set(A, i, 2, c[i]);
    7273  }
    7374  m13 = DetGet(A, 1);
    7475
    7576  for(int i=0;i<3;i++) {
    76     gsl_matrix_set(A, i, 0, a.x[i]*a.x[i] + b.x[i]*b.x[i] + c.x[i]*c.x[i]);
    77     gsl_matrix_set(A, i, 1, a.x[i]);
    78     gsl_matrix_set(A, i, 2, b.x[i]);
     77    gsl_matrix_set(A, i, 0, a[i]*a[i] + b[i]*b[i] + c[i]*c[i]);
     78    gsl_matrix_set(A, i, 1, a[i]);
     79    gsl_matrix_set(A, i, 2, b[i]);
    7980  }
    8081  m14 = DetGet(A, 1);
    … …  
    8384    DoeLog(1) && (eLog()<< Verbose(1) << "three points are colinear." << endl);
    8485
    85   center->x[0] =  0.5 * m12/ m11;
    86   center->x[1] = -0.5 * m13/ m11;
    87   center->x[2] =  0.5 * m14/ m11;
    88 
    89   if (fabs(a.Distance(center) - RADIUS) > MYEPSILON)
    90     DoeLog(1) && (eLog()<< Verbose(1) << "The given center is further way by " << fabs(a.Distance(center) - RADIUS) << " from a than RADIUS." << endl);
     86  center->at(0) =  0.5 * m12/ m11;
     87  center->at(1) = -0.5 * m13/ m11;
     88  center->at(2) =  0.5 * m14/ m11;
     89
     90  if (fabs(a.Distance(*center) - RADIUS) > MYEPSILON)
     91    DoeLog(1) && (eLog()<< Verbose(1) << "The given center is further way by " << fabs(a.Distance(*center) - RADIUS) << " from a than RADIUS." << endl);
    9192
    9293  gsl_matrix_free(A);
    … …  
    120121  Vector OtherCenter;
    121122  Center->Zero();
    122   helper.CopyVector(&a);
    123   helper.Scale(sin(2.*alpha));
    124   Center->AddVector(&helper);
    125   helper.CopyVector(&b);
    126   helper.Scale(sin(2.*beta));
    127   Center->AddVector(&helper);
    128   helper.CopyVector(&c);
    129   helper.Scale(sin(2.*gamma));
    130   Center->AddVector(&helper);
     123  helper = sin(2.*alpha) * a;
     124  (*Center) += helper;
     125  helper = sin(2.*beta) * b;
     126  (*Center) += helper;
     127  helper = sin(2.*gamma) * c;
     128  (*Center) += helper;
    131129  //*Center = a * sin(2.*alpha) + b * sin(2.*beta) + c * sin(2.*gamma) ;
    132130  Center->Scale(1./(sin(2.*alpha) + sin(2.*beta) + sin(2.*gamma)));
    133   NewUmkreismittelpunkt->CopyVector(Center);
     131  (*NewUmkreismittelpunkt) = (*Center);
    134132  DoLog(1) && (Log() << Verbose(1) << "Center of new circumference is " << *NewUmkreismittelpunkt << ".\n");
    135133  // Here we calculated center of circumscribing circle, using barycentric coordinates
    136134  DoLog(1) && (Log() << Verbose(1) << "Center of circumference is " << *Center << " in direction " << *Direction << ".\n");
    137135
    138   TempNormal.CopyVector(&a);
    139   TempNormal.SubtractVector(&b);
    140   helper.CopyVector(&a);
    141   helper.SubtractVector(&c);
    142   TempNormal.VectorProduct(&helper);
     136  TempNormal = a - b;
     137  helper = a - c;
     138  TempNormal.VectorProduct(helper);
    143139  if (fabs(HalfplaneIndicator) < MYEPSILON)
    144140    {
    145       if ((TempNormal.ScalarProduct(AlternativeDirection) <0 && AlternativeIndicator >0) || (TempNormal.ScalarProduct(AlternativeDirection) >0 && AlternativeIndicator <0))
     141      if ((TempNormal.ScalarProduct(*AlternativeDirection) <0 && AlternativeIndicator >0) || (TempNormal.ScalarProduct(*AlternativeDirection) >0 && AlternativeIndicator <0))
    146142        {
    147           TempNormal.Scale(-1);
     143          TempNormal *= -1;
    148144        }
    149145    }
    150146  else
    151147    {
    152       if (((TempNormal.ScalarProduct(Direction)<0) && (HalfplaneIndicator >0)) || ((TempNormal.ScalarProduct(Direction)>0) && (HalfplaneIndicator<0)))
     148      if (((TempNormal.ScalarProduct(*Direction)<0) && (HalfplaneIndicator >0)) || ((TempNormal.ScalarProduct(*Direction)>0) && (HalfplaneIndicator<0)))
    153149        {
    154           TempNormal.Scale(-1);
     150          TempNormal *= -1;
    155151        }
    156152    }
    … …  
    161157  TempNormal.Scale(Restradius);
    162158  DoLog(1) && (Log() << Verbose(1) << "Shift vector to sphere of circumference is " << TempNormal << ".\n");
    163 
    164   Center->AddVector(&TempNormal);
     159  (*Center) += TempNormal;
    165160  DoLog(1) && (Log() << Verbose(1) << "Center of sphere of circumference is " << *Center << ".\n");
    166161  GetSphere(&OtherCenter, a, b, c, RADIUS);
    … …  
    180175  Vector helper;
    181176  double alpha, beta, gamma;
    182   Vector SideA, SideB, SideC;
    183   SideA.CopyVector(b);
    184   SideA.SubtractVector(&c);
    185   SideB.CopyVector(c);
    186   SideB.SubtractVector(&a);
    187   SideC.CopyVector(a);
    188   SideC.SubtractVector(&b);
    189   alpha = M_PI - SideB.Angle(&SideC);
    190   beta = M_PI - SideC.Angle(&SideA);
    191   gamma = M_PI - SideA.Angle(&SideB);
     177  Vector SideA = b - c;
     178  Vector SideB = c - a;
     179  Vector SideC = a - b;
     180  alpha = M_PI - SideB.Angle(SideC);
     181  beta = M_PI - SideC.Angle(SideA);
     182  gamma = M_PI - SideA.Angle(SideB);
    192183  //Log() << Verbose(1) << "INFO: alpha = " << alpha/M_PI*180. << ", beta = " << beta/M_PI*180. << ", gamma = " << gamma/M_PI*180. << "." << endl;
    193184  if (fabs(M_PI - alpha - beta - gamma) > HULLEPSILON) {
    … …  
    196187
    197188  Center->Zero();
    198   helper.CopyVector(a);
    199   helper.Scale(sin(2.*alpha));
    200   Center->AddVector(&helper);
    201   helper.CopyVector(b);
    202   helper.Scale(sin(2.*beta));
    203   Center->AddVector(&helper);
    204   helper.CopyVector(c);
    205   helper.Scale(sin(2.*gamma));
    206   Center->AddVector(&helper);
     189  helper = sin(2.*alpha) * a;
     190  (*Center) += helper;
     191  helper = sin(2.*beta) * b;
     192  (*Center) += helper;
     193  helper = sin(2.*gamma) * c;
     194  (*Center) += helper;
    207195  Center->Scale(1./(sin(2.*alpha) + sin(2.*beta) + sin(2.*gamma)));
    208196};
    … …  
    226214  Vector helper;
    227215  double radius, alpha;
    228   Vector RelativeOldSphereCenter;
    229   Vector RelativeNewSphereCenter;
    230 
    231   RelativeOldSphereCenter.CopyVector(&OldSphereCenter);
    232   RelativeOldSphereCenter.SubtractVector(&CircleCenter);
    233   RelativeNewSphereCenter.CopyVector(&NewSphereCenter);
    234   RelativeNewSphereCenter.SubtractVector(&CircleCenter);
    235   helper.CopyVector(&RelativeNewSphereCenter);
     216
     217  Vector RelativeOldSphereCenter = OldSphereCenter - CircleCenter;
     218  Vector RelativeNewSphereCenter = NewSphereCenter - CircleCenter;
     219  helper = RelativeNewSphereCenter;
    236220  // test whether new center is on the parameter circle's plane
    237   if (fabs(helper.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {
    238     DoeLog(1) && (eLog()<< Verbose(1) << "Something's very wrong here: NewSphereCenter is not on the band's plane as desired by " <<fabs(helper.ScalarProduct(&CirclePlaneNormal))  << "!" << endl);
    239     helper.ProjectOntoPlane(&CirclePlaneNormal);
     221  if (fabs(helper.ScalarProduct(CirclePlaneNormal)) > HULLEPSILON) {
     222    DoeLog(1) && (eLog()<< Verbose(1) << "Something's very wrong here: NewSphereCenter is not on the band's plane as desired by " <<fabs(helper.ScalarProduct(CirclePlaneNormal))  << "!" << endl);
     223    helper.ProjectOntoPlane(CirclePlaneNormal);
    240224  }
    241225  radius = helper.NormSquared();
    … …  
    243227  if (fabs(radius - CircleRadius) > HULLEPSILON)
    244228    DoeLog(1) && (eLog()<< Verbose(1) << "The projected center of the new sphere has radius " << radius << " instead of " << CircleRadius << "." << endl);
    245   alpha = helper.Angle(&RelativeOldSphereCenter);
     229  alpha = helper.Angle(RelativeOldSphereCenter);
    246230  // make the angle unique by checking the halfplanes/search direction
    247   if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)  // acos is not unique on [0, 2.*M_PI), hence extra check to decide between two half intervals
     231  if (helper.ScalarProduct(SearchDirection) < -HULLEPSILON)  // acos is not unique on [0, 2.*M_PI), hence extra check to decide between two half intervals
    248232    alpha = 2.*M_PI - alpha;
    249233  DoLog(1) && (Log() << Verbose(1) << "INFO: RelativeNewSphereCenter is " << helper << ", RelativeOldSphereCenter is " << RelativeOldSphereCenter << " and resulting angle is " << alpha << "." << endl);
    250   radius = helper.Distance(&RelativeOldSphereCenter);
    251   helper.ProjectOntoPlane(&NormalVector);
     234  radius = helper.Distance(RelativeOldSphereCenter);
     235  helper.ProjectOntoPlane(NormalVector);
    252236  // check whether new center is somewhat away or at least right over the current baseline to prevent intersecting triangles
    253237  if ((radius > HULLEPSILON) || (helper.Norm() < HULLEPSILON)) {
    … …  
    279263  struct Intersection *I = (struct Intersection *)params;
    280264  Vector intersection;
    281   Vector SideA,SideB,HeightA, HeightB;
    282265  for (int i=0;i<NDIM;i++)
    283     intersection.x[i] = gsl_vector_get(x, i);
    284 
    285   SideA.CopyVector(&(I->x1));
    286   SideA.SubtractVector(&I->x2);
    287   HeightA.CopyVector(&intersection);
    288   HeightA.SubtractVector(&I->x1);
    289   HeightA.ProjectOntoPlane(&SideA);
    290 
    291   SideB.CopyVector(&I->x3);
    292   SideB.SubtractVector(&I->x4);
    293   HeightB.CopyVector(&intersection);
    294   HeightB.SubtractVector(&I->x3);
    295   HeightB.ProjectOntoPlane(&SideB);
    296 
    297   retval = HeightA.ScalarProduct(&HeightA) + HeightB.ScalarProduct(&HeightB);
     266    intersection[i] = gsl_vector_get(x, i);
     267
     268  Vector SideA = I->x1 -I->x2 ;
     269  Vector HeightA = intersection - I->x1;
     270  HeightA.ProjectOntoPlane(SideA);
     271
     272  Vector SideB = I->x3 - I->x4;
     273  Vector HeightB = intersection - I->x3;
     274  HeightB.ProjectOntoPlane(SideB);
     275
     276  retval = HeightA.ScalarProduct(HeightA) + HeightB.ScalarProduct(HeightB);
    298277  //Log() << Verbose(1) << "MinIntersectDistance called, result: " << retval << endl;
    299278
    … …  
    320299
    321300  struct Intersection par;
    322     par.x1.CopyVector(&point1);
    323     par.x2.CopyVector(&point2);
    324     par.x3.CopyVector(&point3);
    325     par.x4.CopyVector(&point4);
     301    par.x1 = point1;
     302    par.x2 = point2;
     303    par.x3 = point3;
     304    par.x4 = point4;
    326305
    327306    const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex;
    … …  
    336315    /* Starting point */
    337316    x = gsl_vector_alloc(NDIM);
    338     gsl_vector_set(x, 0, point1.x[0]);
    339     gsl_vector_set(x, 1, point1.x[1]);
    340     gsl_vector_set(x, 2, point1.x[2]);
     317    gsl_vector_set(x, 0, point1[0]);
     318    gsl_vector_set(x, 1, point1[1]);
     319    gsl_vector_set(x, 2, point1[2]);
    341320
    342321    /* Set initial step sizes to 1 */
    … …  
    369348
    370349    // check whether intersection is in between or not
    371   Vector intersection, SideA, SideB, HeightA, HeightB;
     350  Vector intersection;
    372351  double t1, t2;
    373352  for (int i = 0; i < NDIM; i++) {
    374     intersection.x[i] = gsl_vector_get(s->x, i);
    375   }
    376 
    377   SideA.CopyVector(&par.x2);
    378   SideA.SubtractVector(&par.x1);
    379   HeightA.CopyVector(&intersection);
    380   HeightA.SubtractVector(&par.x1);
    381 
    382   t1 = HeightA.ScalarProduct(&SideA)/SideA.ScalarProduct(&SideA);
    383 
    384   SideB.CopyVector(&par.x4);
    385   SideB.SubtractVector(&par.x3);
    386   HeightB.CopyVector(&intersection);
    387   HeightB.SubtractVector(&par.x3);
    388 
    389   t2 = HeightB.ScalarProduct(&SideB)/SideB.ScalarProduct(&SideB);
     353    intersection[i] = gsl_vector_get(s->x, i);
     354  }
     355
     356  Vector SideA = par.x2 - par.x1;
     357  Vector HeightA = intersection - par.x1;
     358
     359  t1 = HeightA.ScalarProduct(SideA)/SideA.ScalarProduct(SideA);
     360
     361  Vector SideB = par.x4 - par.x3;
     362  Vector HeightB = intersection - par.x3;
     363
     364  t2 = HeightB.ScalarProduct(SideB)/SideB.ScalarProduct(SideB);
    390365
    391366  Log() << Verbose(1) << "Intersection " << intersection << " is at "
    … …  
    427402  if (point.IsZero())
    428403    return M_PI;
    429   double phi = point.Angle(&reference);
    430   if (OrthogonalVector.ScalarProduct(&point) > 0) {
     404  double phi = point.Angle(reference);
     405  if (OrthogonalVector.ScalarProduct(point) > 0) {
    431406    phi = 2.*M_PI - phi;
    432407  }
    … …  
    451426  double volume;
    452427
    453   TetraederVector[0].CopyVector(a);
    454   TetraederVector[1].CopyVector(b);
    455   TetraederVector[2].CopyVector(c);
     428  TetraederVector[0] = a;
     429  TetraederVector[1] = b;
     430  TetraederVector[2] = c;
    456431  for (int j=0;j<3;j++)
    457     TetraederVector[j].SubtractVector(&d);
    458   Point.CopyVector(&TetraederVector[0]);
    459   Point.VectorProduct(&TetraederVector[1]);
    460   volume = 1./6. * fabs(Point.ScalarProduct(&TetraederVector[2]));
     432    TetraederVector[j].SubtractVector(d);
     433  Point = TetraederVector[0];
     434  Point.VectorProduct(TetraederVector[1]);
     435  volume = 1./6. * fabs(Point.ScalarProduct(TetraederVector[2]));
    461436  return volume;
    462437};
    … …  
    582557        if (List != NULL) {
    583558          for (LinkedCell::LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
    584             helper.CopyVector(Point);
    585             helper.SubtractVector((*Runner)->node);
     559            helper = (*Point) - (*(*Runner)->node);
    586560            double currentNorm = helper. Norm();
    587561            if (currentNorm < distance) {
    … …  
    637611        if (List != NULL) {
    638612          for (LinkedCell::LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
    639             helper.CopyVector(Point);
    640             helper.SubtractVector((*Runner)->node);
     613            helper = (*Point) - (*(*Runner)->node);
    641614            double currentNorm = helper.NormSquared();
    642615            if (currentNorm < distance) {
    … …  
    677650        Info FunctionInfo(__func__);
    678651  // construct the plane of the two baselines (i.e. take both their directional vectors)
    679   Vector Normal;
    680   Vector Baseline, OtherBaseline;
    681   Baseline.CopyVector(Base->endpoints[1]->node->node);
    682   Baseline.SubtractVector(Base->endpoints[0]->node->node);
    683   OtherBaseline.CopyVector(OtherBase->endpoints[1]->node->node);
    684   OtherBaseline.SubtractVector(OtherBase->endpoints[0]->node->node);
    685   Normal.CopyVector(&Baseline);
    686   Normal.VectorProduct(&OtherBaseline);
     652  Vector Baseline = (*Base->endpoints[1]->node->node) - (*Base->endpoints[0]->node->node);
     653  Vector OtherBaseline = (*OtherBase->endpoints[1]->node->node) - (*OtherBase->endpoints[0]->node->node);
     654  Vector Normal = Baseline;
     655  Normal.VectorProduct(OtherBaseline);
    687656  Normal.Normalize();
    688657  DoLog(1) && (Log() << Verbose(1) << "First direction is " << Baseline << ", second direction is " << OtherBaseline << ", normal of intersection plane is " << Normal << "." << endl);
    689658
    690659  // project one offset point of OtherBase onto this plane (and add plane offset vector)
    691   Vector NewOffset;
    692   NewOffset.CopyVector(OtherBase->endpoints[0]->node->node);
    693   NewOffset.SubtractVector(Base->endpoints[0]->node->node);
    694   NewOffset.ProjectOntoPlane(&Normal);
    695   NewOffset.AddVector(Base->endpoints[0]->node->node);
    696   Vector NewDirection;
    697   NewDirection.CopyVector(&NewOffset);
    698   NewDirection.AddVector(&OtherBaseline);
     660  Vector NewOffset = (*OtherBase->endpoints[0]->node->node) - (*Base->endpoints[0]->node->node);
     661  NewOffset.ProjectOntoPlane(Normal);
     662  NewOffset += (*Base->endpoints[0]->node->node);
     663  Vector NewDirection = NewOffset + OtherBaseline;
    699664
    700665  // calculate the intersection between this projected baseline and Base
    701666  Vector *Intersection = new Vector;
    702   Intersection->GetIntersectionOfTwoLinesOnPlane(Base->endpoints[0]->node->node, Base->endpoints[1]->node->node, &NewOffset, &NewDirection, &Normal);
    703   Normal.CopyVector(Intersection);
    704   Normal.SubtractVector(Base->endpoints[0]->node->node);
    705   DoLog(1) && (Log() << Verbose(1) << "Found closest point on " << *Base << " at " << *Intersection << ", factor in line is " << fabs(Normal.ScalarProduct(&Baseline)/Baseline.NormSquared()) << "." << endl);
     667  *Intersection = GetIntersectionOfTwoLinesOnPlane(*(Base->endpoints[0]->node->node),
     668                                                   *(Base->endpoints[1]->node->node),
     669                                                     NewOffset, NewDirection);
     670  Normal = (*Intersection) - (*Base->endpoints[0]->node->node);
     671  DoLog(1) && (Log() << Verbose(1) << "Found closest point on " << *Base << " at " << *Intersection << ", factor in line is " << fabs(Normal.ScalarProduct(Baseline)/Baseline.NormSquared()) << "." << endl);
    706672
    707673  return Intersection;
    … …  
    721687    return -1;
    722688  }
    723   distance = x->DistanceToPlane(&triangle->NormalVector, triangle->endpoints[0]->node->node);
     689  distance = x->DistanceToPlane(triangle->NormalVector, *triangle->endpoints[0]->node->node);
    724690  return distance;
    725691};
    … …  
    747713      *vrmlfile << "Sphere {" << endl << "  "; // 2 is sphere type
    748714      for (i=0;i<NDIM;i++)
    749         *vrmlfile << Walker->node->x[i]-center->x[i] << " ";
     715        *vrmlfile << Walker->node->at(i)-center->at(i) << " ";
    750716      *vrmlfile << "\t0.1\t1. 1. 1." << endl; // radius 0.05 and white as colour
    751717      cloud->GoToNext();
    … …  
    757723      for (i=0;i<3;i++) { // print each node
    758724        for (int j=0;j<NDIM;j++)  // and for each node all NDIM coordinates
    759           *vrmlfile << TriangleRunner->second->endpoints[i]->node->node->x[j]-center->x[j] << " ";
     725          *vrmlfile << TriangleRunner->second->endpoints[i]->node->node->at(j)-center->at(j) << " ";
    760726        *vrmlfile << "\t";
    761727      }
    … …  
    784750    Vector *center = cloud->GetCenter();
    785751    // make the circumsphere's center absolute again
    786     helper.CopyVector(Tess->LastTriangle->endpoints[0]->node->node);
    787     helper.AddVector(Tess->LastTriangle->endpoints[1]->node->node);
    788     helper.AddVector(Tess->LastTriangle->endpoints[2]->node->node);
    789     helper.Scale(1./3.);
    790     helper.SubtractVector(center);
     752    Vector helper = (1./3.) * ((*Tess->LastTriangle->endpoints[0]->node->node) +
     753                               (*Tess->LastTriangle->endpoints[1]->node->node) +
     754                               (*Tess->LastTriangle->endpoints[2]->node->node));
     755    helper -= (*center);
    791756    // and add to file plus translucency object
    792757    *rasterfile << "# current virtual sphere\n";
    793758    *rasterfile << "8\n  25.0    0.6     -1.0 -1.0 -1.0     0.2        0 0 0 0\n";
    794     *rasterfile << "2\n  " << helper.x[0] << " " << helper.x[1] << " " << helper.x[2] << "\t" << 5. << "\t1 0 0\n";
     759    *rasterfile << "2\n  " << helper[0] << " " << helper[1] << " " << helper[2] << "\t" << 5. << "\t1 0 0\n";
    795760    *rasterfile << "9\n  terminating special property\n";
    796761    delete(center);
    … …  
    820785      *rasterfile << "2" << endl << "  ";  // 2 is sphere type
    821786      for (i=0;i<NDIM;i++)
    822         *rasterfile << Walker->node->x[i]-center->x[i] << " ";
     787        *rasterfile << Walker->node->at(i)-center->at(i) << " ";
    823788      *rasterfile << "\t0.1\t1. 1. 1." << endl; // radius 0.05 and white as colour
    824789      cloud->GoToNext();
    … …  
    831796      for (i=0;i<3;i++) {  // print each node
    832797        for (int j=0;j<NDIM;j++)  // and for each node all NDIM coordinates
    833           *rasterfile << TriangleRunner->second->endpoints[i]->node->node->x[j]-center->x[j] << " ";
     798          *rasterfile << TriangleRunner->second->endpoints[i]->node->node->at(j)-center->at(j) << " ";
    834799        *rasterfile << "\t";
    835800      }
    … …  
    881846      Walker = target->second->node;
    882847      LookupList[Walker->nr] = Counter++;
    883       *tecplot << Walker->node->x[0] << " " << Walker->node->x[1] << " " << Walker->node->x[2] << " " << target->second->value << endl;
     848      *tecplot << Walker->node->at(0) << " " << Walker->node->at(1) << " " << Walker->node->at(2) << " " << target->second->value << endl;
    884849    }
    885850    *tecplot << endl;
Note: See TracChangeset for help on using the changeset viewer.