Added a space tweezer with report of position, frequency and amplitude - to be fully tested

.gitignore

@@ -18,6 +18,9 @@ bin3
 *.mode1v3
 DerivedData
 
+# Xcode
+cytosim.xcodeproj/xcuserdata/
+cytosim.xcodeproj/xcshareddata/xcschemes/
 
 # Created by https://www.gitignore.io/api/c++,cmake
 

src/sim/makefile.inc

@@ -6,7 +6,7 @@ OBJ_SPACE := space.o space_prop.o space_square.o space_sphere.o\
              space_dice.o space_torus.o space_polygon.o\
              space_polygonZ.o space_ellipse.o space_cylinder.o space_ring.o\
              space_cylinderZ.o space_capsule.o space_strip.o space_periodic.o\
-             space_banana.o space_cylinderP.o
+             space_banana.o space_cylinderP.o space_tweezer.o
 
 OBJ_HANDS := hand.o hand_prop.o\
              motor.o motor_prop.o\

src/sim/simul.h

@@ -517,6 +517,9 @@ public:
     /// a special print for Romain Gibeaux
     void      reportAshbya(std::ostream&) const;
     
+    /// print tweezer position
+    void      reportTweezer(std::ostream&, std::string const&) const;
+    
     /// print something
     void      reportCustom(std::ostream&) const;
 

src/sim/simul_report.cc

@@ -161,6 +161,7 @@ void Simul::report(std::ostream& out, std::string arg, Glossary& opt) const
  `couple:anatomy`        | Composition of couples
  `couple:NAME`           | Position of couples of class NAME
  `couple:hands`          | Composition of couples
+ `space:tweezer`     | Position of the tweezer
  
  */
 void Simul::report0(std::ostream& out, std::string const& arg, Glossary& opt) const
@@ -263,6 +264,8 @@ void Simul::report0(std::ostream& out, std::string const& arg, Glossary& opt) co
     {
         if ( what == "force" )
             return reportSpaceForce(out);
+        else if ( what == "tweezer" )
+            return reportTweezer(out, what);
         else if ( what.empty() )
             return reportSpace(out);
         throw InvalidSyntax("I only know `space'");
@@ -1391,6 +1394,29 @@ void Simul::reportSpaceForce(std::ostream& out) const
 }
 
 
+/**
+ Export the position of a  tweezer Space
+ */
+
+void Simul::reportTweezer(std::ostream& out, std::string const& what) const
+{
+    if ( properties.find("space", what) == nullptr )
+        throw InvalidParameter("Unknown space `"+what+"'");
+
+    for ( Space * obj=spaces.first(); obj; obj=obj->next() )
+    {
+        if ( obj->property()->name() == what )
+        {
+        out << LIN << obj->prop->name();
+        out << SEP << obj->identity();
+        out << SEP << obj->prop->shape;
+        out << '\n' ;
+        obj->report(out);
+        }
+    }
+}
+
+
 /**
  Report quantity of substance in Field
  */

src/sim/space.cc

@@ -391,6 +391,12 @@ void Space::write(Outputter& out) const
     out.writeUInt16(0);
 }
 
+void Space::report(std::ostream & out) const
+{
+    //std::cerr << "Space : " << prop->shape << std::endl;
+    throw InvalidIO("Space::report unimplemented");
+}
+
 
 void Space::read(Inputter& in, Simul&, ObjectTag)
 {

src/sim/space.h

@@ -176,6 +176,9 @@ public:
     /// read from file
     void           read(Inputter&, Simul&, ObjectTag);
     
+    /// report to file (unimplemented by default)
+    virtual void  report(std::ostream &) const;
+    
     /// get dimensions from array `len`
     virtual void   setLengths(const real len[8]) {}
 

src/sim/space_prop.cc

@@ -23,6 +23,7 @@
 #include "space_cylinderZ.h"
 #include "space_cylinderP.h"
 #include "space_ring.h"
+#include "space_tweezer.h"
 
 /**
  @defgroup SpaceGroup Space and Geometry
@@ -64,6 +65,7 @@
  `cylinderZ`   | SpaceCylinderZ       | radius bottom top
  `cylinderP`   | SpaceCylinderP       | length radius
  `ring`        | SpaceRing            | length radius
+ `tweezer`     | SpaceTweezer     | centerX centerY centerZ amplitude frequency
  
  Example:
  
@@ -92,6 +94,7 @@ Space * SpaceProp::newSpace() const
     if ( s=="strip" || s=="semi_periodic" )    return new SpaceStrip(this);
     if ( s=="periodic" )                       return new SpacePeriodic(this);
     if ( s=="ellipse" || s=="ellipsoid" )      return new SpaceEllipse(this);
+    if ( s=="tweezer" )                        return new SpaceTweezer(this);
 #if ( DIM >= 3 )
     if ( s=="cubic" )                          return new SpaceSquare(this);
     if ( s=="cylinder" )                       return new SpaceCylinder(this);
@@ -180,7 +183,6 @@ void SpaceProp::read(Glossary& glos)
         }
 #endif
     }
-    
     if ( glos.set(display, "display") )
         display_fresh = true;
 }
@@ -191,6 +193,9 @@ void SpaceProp::complete(Simul const& sim)
 {
     if ( shape.empty() )
         throw InvalidParameter("space:shape must be defined");
+    
+    // initializes the time step here
+    time_step = sim.time_step();
 }
 
 //------------------------------------------------------------------------------
@@ -200,5 +205,7 @@ void SpaceProp::write_values(std::ostream& os) const
     //write_value(os, "geometry",   geometry);
     write_value(os, "shape",      shape);
     write_value(os, "display",    "("+display+")");
+//    if ( tweezer_space )
+//        write_value(os, "frequency",  frequency);
 }
 

src/sim/space_prop.h

@@ -39,6 +39,9 @@ public:
     
     /// BACKWARD_COMPATIBILITY
     std::string  dimensions_;
+        
+    /// time step for dynamic spaces
+    real time_step;
     
     /// derived variable: flag to indicate that `display` has a new value
     bool         display_fresh;

src/sim/spaces/CMakeLists.txt

@@ -15,9 +15,9 @@ set(SPACES_SOURCES
     "${PROJECT_SOURCE_DIR}/src/sim/spaces/space_strip.cc"
     "${PROJECT_SOURCE_DIR}/src/sim/spaces/space_periodic.cc"
     "${PROJECT_SOURCE_DIR}/src/sim/spaces/space_cylinderP.cc"
-	
     "${PROJECT_SOURCE_DIR}/src/sim/spaces/space_polygonZ.cc"
     "${PROJECT_SOURCE_DIR}/src/sim/spaces/space_ring.cc"
+    "${PROJECT_SOURCE_DIR}/src/sim/spaces/space_tweezer.cc"
 )
 
 add_library(${SPACES_LIB_TARGET} STATIC ${SPACES_SOURCES})

src/sim/spaces/space_tweezer.cc

+// Cytosim 3.0
+//
+//  File created by Hervé Turlier on 2/13/19.
+//  Modified on 1/14/22
+//  Copyright © 2019 Collège-de-France. All rights reserved.
+
+#include "space_tweezer.h"
+#include "exceptions.h"
+#include "smath.h"
+#include "simul.h"
+#include <math.h>
+#include "meca.h"
+#include "iowrapper.h"
+
+
+SpaceTweezer::SpaceTweezer(const SpaceProp* p) :
+Space(p)
+{
+    amplitude_ = 0;
+    center_.reset();
+    frequency_ = 0;
+    time = 0;
+}
+
+//------------------------------------------------------------------------------
+//                         Geometrical redefinitions
+//------------------------------------------------------------------------------
+
+void SpaceTweezer::resize(Glossary& opt)
+{
+    real amp = amplitude_;
+    real freq = frequency_;
+    Vector cent = center_;
+    
+//    for ( int d = 0; d < DIM; ++d )
+//        opt.set(center_[d], "center", d);
+
+    opt.set(cent, "center");
+    center_ = cent;
+
+    opt.set(amp, "amplitude");
+    if ( amp < 0 )
+        throw InvalidParameter(prop->name()+"amplitude should be >= 0");
+    amplitude_ = amp;
+
+    opt.set(freq, "frequency");
+    if ( freq < 0 )
+        throw InvalidParameter(prop->name()+"frequency should be >= 0");
+    frequency_ = freq;
+        
+    update();
+}
+
+void SpaceTweezer::update()
+{
+    // nothing to do
+}
+
+void SpaceTweezer::boundaries(Vector& inf, Vector& sup) const
+{
+    inf.set(-amplitude_, 0, 0);
+    sup.set(amplitude_, 0, 0);
+}
+
+
+real SpaceTweezer::volume() const
+{
+    return 0;
+}
+
+// return false to always have an interaction
+bool SpaceTweezer::inside(Vector const& pos) const
+{
+    return false;
+}
+
+// return false to always have an interaction
+bool  SpaceTweezer::allInside( Vector const&, real rad ) const
+{
+    return false;
+}
+
+// return true to always have an interaction
+bool  SpaceTweezer::allOutside( Vector const&, real rad ) const
+{
+    return true;
+}
+
+Vector  SpaceTweezer::project(Vector const& pos) const
+{
+    return center_;
+}
+
+//------------------------------------------------------------------------------
+//                         Interaction -> trap force
+//------------------------------------------------------------------------------
+
+
+void SpaceTweezer::setInteraction(Vector const& pos, Mecapoint const& pe, Meca & meca, real stiff) const
+{
+    meca.addPointClamp(pe, center_, stiff);
+}
+
+
+void SpaceTweezer::setInteraction(Vector const& pos, Mecapoint const& pe, real rad, Meca & meca, real stiff) const
+{
+    meca.addPointClamp(pe, center_, stiff);
+}
+
+
+//------------------------------------------------------------------------------
+//                                  IO
+//------------------------------------------------------------------------------
+
+void  SpaceTweezer::write(Outputter& out) const
+{
+    out.put_characters("tweezer", 16);
+    out.writeUInt16(2+DIM);
+    out.writeFloat(amplitude_);
+    out.writeFloat(frequency_);
+    for ( int d = 0; d < DIM ; ++d )
+        out.writeDouble(center_[d]);
+}
+
+void SpaceTweezer::setLengths(const real len[])
+{
+    amplitude_ = len[0];
+    frequency_ = len[1];
+    for (int d = 0; d < DIM ; ++d)
+        center_[d] = len[2+d];
+    update();
+}
+
+void  SpaceTweezer::read(Inputter& in, Simul& sim, ObjectTag)
+{
+    real len[8] = { 0 };
+    read_data(in, len, "tweezer");
+    setLengths(len);
+}
+
+void SpaceTweezer::report(std::ostream& out) const
+{
+    out << "% tweezer position" << '\n';
+    out << center_ ;
+    out << '\n';
+    out << "% tweezer frequency" << '\n';
+    out << frequency_ ;
+    out << '\n';
+    out << "% tweezer amplitude" << '\n';
+    out << amplitude_ ;
+    out << '\n';
+}
+
+//------------------------------------------------------------------------------
+//                      Dynamic space time-stepping
+//------------------------------------------------------------------------------
+
+
+void SpaceTweezer::step()
+{
+    // updating time
+    time += prop->time_step;
+    // updating trap position
+    center_.XX = amplitude_*cos(2*M_PI*frequency_*time + M_PI/2);
+}
+
+
+//------------------------------------------------------------------------------
+//                         OPENGL  DISPLAY
+//------------------------------------------------------------------------------
+
+#ifdef DISPLAY
+#include "opengl.h"
+#include "gle.h"
+
+
+bool SpaceTweezer::draw() const
+{
+#if ( DIM > 2 )
+
+    const size_t fin = 512;
+
+    GLfloat L = (GLfloat)amplitude_;
+    GLfloat R = (GLfloat)0.1;
+    
+    GLfloat c[fin+1], s[fin+1];
+    gle::circle(fin, c, s, 1);
+    
+    glBegin(GL_TRIANGLE_STRIP);
+    for ( size_t sc = 0; sc <= fin; ++sc )
+    {
+        GLfloat ca = c[sc], sa = s[sc];
+        glNormal3f( 0, ca, sa );
+        glVertex3f( +L, R*ca, R*sa );
+        glVertex3f( -L, R*ca, R*sa );
+    }
+    glEnd();
+    
+    // draw the cap:
+    glBegin(GL_TRIANGLE_FAN);
+    glNormal3f( +1, 0, 0 );
+    glVertex3f( +L, 0, 0 );
+    for ( size_t sc = 0; sc <= fin; ++sc )
+        glVertex3f( +L, R*c[sc], R*s[sc] );
+    glEnd();
+    
+    // draw the cap:
+    glBegin(GL_TRIANGLE_FAN);
+    glNormal3f( -1, 0, 0 );
+    glVertex3f( -L, 0, 0 );
+    for ( size_t sc = 0; sc <= fin; ++sc )
+        glVertex3f( -L,-R*c[sc], R*s[sc] );
+    glEnd();
+    
+#endif
+    return true;
+}
+
+#else
+
+bool SpaceTweezer::draw() const
+{
+    return false;
+}
+
+#endif

src/sim/spaces/space_tweezer.h

+// Cytosim 3.0
+//
+//  File created by Hervé Turlier on 2/13/19.
+//  Copyright © 2019 Collège-de-France. All rights reserved.
+
+#ifndef space_tweezer_h
+#define space_tweezer_h
+
+#include "dim.h"
+#include "space_sphere.h"
+
+/// An harmonic trap under sinusoidal movement to simulate an (optical) tweezer
+/**
+ Space `tweezer` is a point-like space (zero volume), centered around the origin, under sinusoidal movement of given amplitude in x-direction
+ An object confined in this spaces will feel a linear restoring force characterized by the 'stiffness' of the interaction.
+ 
+ - amplitude = amplitude of the oscillatory movement of the trap in x-direction
+ - frequency = driving frequency of the oscillatory trap movement
+
+ @ingroup SpaceGroup
+ Forces registered with 'setInteractions' are added, and used to create the trap force.
+ The harmonic trap stiffness is set directly in the interaction.
+ The trap force on an object can be calculated a posteriori from the interaction stiffness k, trap position x_trap and trapped object center of mass x_CM as
+ @code
+ F = k (x_CM - x_trap)
+ @endcode
+
+ To extract the trap position over time, the report function can be used with argument space:tweezer.
+ @ingroup SpaceGroup
+ 
+ 
+ Created by Hervé Turlier, Paris 13.02.2019
+ Modified on 03.03.2019
+ */
+
+class SpaceTweezer : public Space
+{
+
+private:
+    
+    /// trap center
+    Vector      center_;
+
+    /// frequency of trap harmonic driving
+    real        frequency_;
+    
+    /// the amplitude of trap driving
+    real        amplitude_;
+    //real        amplitude() const { return Space::length(0); }
+    
+    /// simulation time
+    real        time;
+
+public:
+    /// constructor
+    SpaceTweezer(const SpaceProp*);
+    
+    /// check number and validity of specified lengths
+    // void        resize() { Space::checkLengths(1, true); }
+    void        resize(Glossary& opt);
+
+    /// no update implemented at this point
+    void  update();
+    
+    // maximum extension along each axis
+    //Vector      extension() const;
+    
+    /// return bounding box in `inf` and `sup`
+    void        boundaries(Vector& inf, Vector& sup) const;
+    
+    /// the volume inside (zero)
+    real        volume() const;
+    
+    /// true if the point is inside the Space - always true here to have an interaction
+    //bool        inside(const real point[]) const;
+    bool        inside(Vector const&) const;
+
+    /// true if the point is inside the Space - always true here to have an interaction
+    //bool        allInside(const real point[], real rad) const;
+    bool        allInside(Vector const&, real rad) const;
+
+    /// always true here to have an interaction
+    //bool        allOutside(const real point[], real rad) const;
+    bool        allOutside(Vector const&, real rad) const;
+
+    /// set `proj` as the point on the edge that is closest to `point`
+    Vector      project(Vector const& pos) const;
+
+    /// apply an harmonic force directed towards the center of the trap
+    //void        setInteraction(Vector const& pos, PointExact const&, Meca &, real stiff) const;
+    void        setInteraction(Vector const& pos, Mecapoint const&, Meca &, real stiff) const;
+
+    /// apply an harmonic force directed towards the center of the trap
+    //void        setInteraction(Vector const& pos, PointExact const&, real rad, Meca &, real stiff) const;
+    void        setInteraction(Vector const& pos, Mecapoint const&, real rad, Meca &, real stiff) const;
+
+    /// the step function will move the trap center
+    void        step();
+    
+    /// read from file
+    //virtual void  read(InputWrapper& , Simul&);
+    void        read(Inputter&, Simul&, ObjectTag);
+
+    /// write to file
+    //virtual void  write(OutputWrapper&) const;
+    void        write(Outputter&) const;
+
+    /// get dimensions from array 'len'
+    void        setLengths(const real len[8]);
+    
+    /// return tweezer center
+    void report(std::ostream&) const;
+    
+    /// OpenGL display function, return true is display was done
+    //bool        display() const;
+    bool        draw() const;
+
+};
+#endif /* space_tweezer_h */