Changeset 92bc5c for src/Fragmentation


Ignore:
Timestamp:
Sep 12, 2016, 2:03:15 PM (8 years ago)
Author:
Frederik Heber <heber@…>
Branches:
Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_StructOpt_integration_tests, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, 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_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, GeometryObjects, Gui_displays_atomic_force_velocity, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, Ubuntu_1604_changes, stable
Children:
8e4471
Parents:
f9d85f
git-author:
Frederik Heber <heber@…> (06/05/14 17:17:10)
git-committer:
Frederik Heber <heber@…> (09/12/16 14:03:15)
Message:

Adding QuaternionTest for SphericalPointDistributionUnitTest, fixed tests.

Location:
src/Fragmentation/Exporters/unittests
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • src/Fragmentation/Exporters/unittests/SphericalPointDistributionUnitTest.cpp

    rf9d85f r92bc5c  
    4646
    4747#include <boost/assign.hpp>
     48#include <boost/math/quaternion.hpp>
    4849
    4950#include "CodePatterns/Assert.hpp"
    5051#include "CodePatterns/Log.hpp"
     52
     53#include "LinearAlgebra/Line.hpp"
    5154
    5255#include "Fragmentation/Exporters/SphericalPointDistribution.hpp"
    … …  
    7780}
    7881
     82void SphericalPointDistributionTest::QuaternionTest()
     83{
     84  Vector oldCenter(0.,1.,0.);
     85  Vector newCenter(1.,0.,0.);
     86
     87  {
     88    // setup quaternion
     89    Vector RotationAxis = newCenter;
     90    RotationAxis.VectorProduct(oldCenter);
     91    RotationAxis.Normalize();
     92    const double RotationAngle = oldCenter.Angle(newCenter)/(M_PI/2.);
     93  //            RotationAxis.Angle(oldCenter) - RotationAxis.Angle(newCenter);
     94    boost::math::quaternion<double> q
     95        (RotationAngle, RotationAxis[0], RotationAxis[1], RotationAxis[2]);
     96    LOG(5, "DEBUG: RotationAxis is " << RotationAxis
     97        << ", RotationAngle is " << RotationAngle);
     98    LOG(5, "DEBUG: Quaternion describing rotation is " << q);
     99    boost::math::quaternion<double> q_inverse =
     100        boost::math::conj(q)/(boost::math::norm(q));
     101    LOG(5, "DEBUG: Quaternion inverse is " << q_inverse);
     102    boost::math::quaternion<double> identity(1,0,0,0);
     103    const boost::math::quaternion<double> unity = q*q_inverse;
     104    LOG(5, "DEBUG: q * q^-1 is " << unity);
     105    CPPUNIT_ASSERT( boost::math::norm(unity - identity) < std::numeric_limits<double>::epsilon()*1e4);
     106
     107    // check that rotation works
     108    boost::math::quaternion<double> p(0., newCenter[0], newCenter[1], newCenter[2]);
     109    LOG(5, "DEBUG: Original newCenter is " << p);
     110    p = p * q_inverse;
     111    p = q * p;
     112    LOG(5, "DEBUG: Rotated newCenter is " << p);
     113    boost::math::quaternion<double> comparison(0., -oldCenter[0], oldCenter[1], oldCenter[2]);
     114    LOG(5, "DEBUG: Difference norm is " << boost::math::norm(p - comparison));
     115    CPPUNIT_ASSERT( boost::math::norm(p - comparison) < std::numeric_limits<double>::epsilon()*1e4);
     116  }
     117
     118  // rotating with angle = 0 flips the vector unwantedly
     119  {
     120    // setup quaternion
     121    Vector RotationAxis = newCenter;
     122    RotationAxis.VectorProduct(oldCenter);
     123    RotationAxis.Normalize();
     124    const double RotationAngle = 0.;
     125  //            RotationAxis.Angle(oldCenter) - RotationAxis.Angle(newCenter);
     126    boost::math::quaternion<double> q
     127        (RotationAngle, RotationAxis[0], RotationAxis[1], RotationAxis[2]);
     128    LOG(5, "DEBUG: RotationAxis is " << RotationAxis
     129        << ", RotationAngle is " << RotationAngle);
     130    LOG(5, "DEBUG: Quaternion describing rotation is " << q);
     131    boost::math::quaternion<double> q_inverse =
     132        boost::math::conj(q)/(boost::math::norm(q));
     133    LOG(5, "DEBUG: Quaternion inverse is " << q_inverse);
     134    boost::math::quaternion<double> identity(1,0,0,0);
     135    const boost::math::quaternion<double> unity = q*q_inverse;
     136    LOG(5, "DEBUG: q * q^-1 is " << unity);
     137    CPPUNIT_ASSERT( boost::math::norm(unity - identity) < std::numeric_limits<double>::epsilon()*1e4);
     138
     139    // check that rotation works
     140    boost::math::quaternion<double> p(0., newCenter[0], newCenter[1], newCenter[2]);
     141    boost::math::quaternion<double> comparison(0., -newCenter[0], newCenter[1], newCenter[2]);
     142    LOG(5, "DEBUG: Original newCenter is " << p);
     143    p = p * q_inverse;
     144    p = q * p;
     145    LOG(5, "DEBUG: Rotated newCenter is " << p);
     146    LOG(5, "DEBUG: Difference norm is " << boost::math::norm(p - comparison));
     147    CPPUNIT_ASSERT( boost::math::norm(p - comparison) < std::numeric_limits<double>::epsilon()*1e4);
     148  }
     149}
     150
    79151/** UnitTest for matchSphericalPointDistributions() with two points
    80152 */
    … …  
    126198    CPPUNIT_ASSERT_EQUAL( expected, remaining );
    127199  }
     200
     201  // test with one point, full rotation
     202  {
     203    Line RotationAxis(zeroVec, Vector(0.2, 0.43, 0.6893248));
     204    SphericalPointDistribution::Polygon_t polygon;
     205    polygon += RotationAxis.rotateVector(Vector(1.,0.,0.), 47.6/180*M_PI);
     206    SphericalPointDistribution::Polygon_t newpolygon =
     207        SPD.get<2>();
     208    SphericalPointDistribution::Polygon_t expected;
     209    expected += RotationAxis.rotateVector(Vector(-1.,0.,0.), 47.6/180*M_PI);
     210    SphericalPointDistribution::Polygon_t remaining =
     211        SphericalPointDistribution::matchSphericalPointDistributions(
     212            polygon,
     213            newpolygon);
     214    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     215  }
    128216}
    129217
    … …  
    149237  }
    150238
    151   // test with one point, just a flip of axis
     239  // test with one point, just a flip of x and y axis
    152240  {
    153241    SphericalPointDistribution::Polygon_t polygon;
    … …  
    157245    SphericalPointDistribution::Polygon_t expected = newpolygon;
    158246    expected.pop_front(); // remove first point
     247    for (SphericalPointDistribution::Polygon_t::iterator iter = expected.begin();
     248        iter != expected.end(); ++iter) {
     249      std::swap((*iter)[0], (*iter)[1]);
     250      (*iter)[0] *= -1.;
     251    }
    159252    SphericalPointDistribution::Polygon_t remaining =
    160253        SphericalPointDistribution::matchSphericalPointDistributions(
    … …  
    194287    SphericalPointDistribution::Polygon_t expected = newpolygon;
    195288    expected.pop_front(); // remove first point
     289    for (SphericalPointDistribution::Polygon_t::iterator iter = expected.begin();
     290        iter != expected.end(); ++iter) {
     291      std::swap((*iter)[0], (*iter)[1]);
     292      (*iter)[0] *= -1.;
     293    }
    196294    SphericalPointDistribution::Polygon_t remaining =
    197295        SphericalPointDistribution::matchSphericalPointDistributions(
    … …  
    231329    SphericalPointDistribution::Polygon_t expected = newpolygon;
    232330    expected.pop_front(); // remove first point
     331    for (SphericalPointDistribution::Polygon_t::iterator iter = expected.begin();
     332        iter != expected.end(); ++iter) {
     333      std::swap((*iter)[0], (*iter)[1]);
     334      (*iter)[0] *= -1.;
     335    }
    233336    SphericalPointDistribution::Polygon_t remaining =
    234337        SphericalPointDistribution::matchSphericalPointDistributions(
    … …  
    268371    SphericalPointDistribution::Polygon_t expected = newpolygon;
    269372    expected.pop_front(); // remove first point
     373    for (SphericalPointDistribution::Polygon_t::iterator iter = expected.begin();
     374        iter != expected.end(); ++iter) {
     375      std::swap((*iter)[0], (*iter)[1]);
     376      (*iter)[0] *= -1.;
     377    }
    270378    SphericalPointDistribution::Polygon_t remaining =
    271379        SphericalPointDistribution::matchSphericalPointDistributions(
    … …  
    305413    SphericalPointDistribution::Polygon_t expected = newpolygon;
    306414    expected.pop_front(); // remove first point
     415    for (SphericalPointDistribution::Polygon_t::iterator iter = expected.begin();
     416        iter != expected.end(); ++iter) {
     417      std::swap((*iter)[0], (*iter)[1]);
     418      (*iter)[0] *= -1.;
     419    }
    307420    SphericalPointDistribution::Polygon_t remaining =
    308421        SphericalPointDistribution::matchSphericalPointDistributions(
    … …  
    342455    SphericalPointDistribution::Polygon_t expected = newpolygon;
    343456    expected.pop_front(); // remove first point
    344     SphericalPointDistribution::Polygon_t remaining =
    345         SphericalPointDistribution::matchSphericalPointDistributions(
    346             polygon,
    347             newpolygon);
    348     CPPUNIT_ASSERT_EQUAL( expected, remaining );
    349   }
    350 }
     457    for (SphericalPointDistribution::Polygon_t::iterator iter = expected.begin();
     458        iter != expected.end(); ++iter) {
     459      std::swap((*iter)[0], (*iter)[1]);
     460      (*iter)[0] *= -1.;
     461    }
     462   SphericalPointDistribution::Polygon_t remaining =
     463        SphericalPointDistribution::matchSphericalPointDistributions(
     464            polygon,
     465            newpolygon);
     466    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     467  }
     468}

  • src/Fragmentation/Exporters/unittests/SphericalPointDistributionUnitTest.hpp

    rf9d85f r92bc5c  
    2222{
    2323    CPPUNIT_TEST_SUITE( SphericalPointDistributionTest) ;
     24    CPPUNIT_TEST ( QuaternionTest );
    2425    CPPUNIT_TEST ( matchSphericalPointDistributionsTest_2 );
    2526    CPPUNIT_TEST ( matchSphericalPointDistributionsTest_3 );
    … …  
    3435      void setUp();
    3536      void tearDown();
     37      void QuaternionTest();
    3638      void matchSphericalPointDistributionsTest_2();
    3739      void matchSphericalPointDistributionsTest_3();
Note: See TracChangeset for help on using the changeset viewer.