Merge branch 'master' into Milan

.gitignore

@@ -27,3 +27,4 @@ cpp/Build-*
 cpp/build
 cpp/build-*
 .aider*
+.env

cpp/common/Interfaces.h

@@ -0,0 +1,64 @@
+#ifndef  Interfaces_h
+#define  Interfaces_h
+
+#include "Vec2.h"
+#include "Vec3.h"
+
+class AnyControler{ public:
+    //virtual void control(void* obj, double dt)=0;
+    virtual void update(double dt)=0;
+};
+
+// Interf
+/**
+ * @brief The BindLoader class is an abstract base class that defines the interface for objects that can bind another object and load data from it
+ */
+class BindLoader {
+public:
+    /**
+     * @brief Binds and loads data from specified object.
+     * @param o A pointer to the object to bind and load data from
+     * @return An integer indicating the result of the operation.
+     */
+    virtual int bindLoad(void* o) = 0;
+};
+
+/**
+ * @brief The Picker class is an abstract base class that defines the interface for picking operations.
+ */
+class Picker {
+public:
+    /**
+     * @brief Picks the nearest object along the given ray.
+     * 
+     * @param ray0 The starting point of the ray.
+     * @param hray The direction and length of the ray.
+     * @param ipick The index of the picked object.
+     * @param mask The mask specifying which objects to consider for picking.
+     * @return The type of the picked object.
+     */
+    virtual int pick_nearest(Vec3d ray0, Vec3d hray, int& ipick, int mask = 0xFFFFFFFF, double Rmax=0.0 ) = 0;
+
+    /**
+     * @brief Picks all objects intersected by the given ray.
+     * 
+     * @param ray0 The starting point of the ray.
+     * @param hray The direction and length of the ray.
+     * @param out  The output array of picked objects. 
+     * @param mask The mask specifying which objects to consider for picking.
+     * @return The number of objects added to the picking list.
+     */
+    virtual int pick_all(Vec3d ray0, Vec3d hray, int* out, int mask = 0xFFFFFFFF, double Rmax=0.0 ) = 0;
+
+    /**
+     * @brief Gets the object associated with the specified index in the picking list.
+     * 
+     * @param picked The index of the picked object.
+     * @param mask   The mask specifying which objects to consider for picking.
+     * @return       A pointer to the picked object.
+     */
+    virtual void* getPickedObject(int picked, int mask = 0xFFFFFFFF ) = 0;
+
+};
+
+#endif

cpp/common/molecular/GridFF.h

@@ -998,6 +998,30 @@ inline void addForce( const Vec3d& pos, const Quat4f& PLQ, Quat4f& fe ) const {
         printf( "GridFF::makeGridFF_Bspline_d() DONE\n" );
     }
 
+
+
+#ifdef WITH_FFTW
+
+    void makeCoulombEwald(  int natoms_, Vec3d * apos_, Quat4d * REQs_, double* VCoul=0 ){
+        std::vector<double> qs(natoms_); 
+        for( int i=0; i<natoms_; i++ ){ qs[i] = REQs_[i].z; }
+        ewald->copy( grid );
+        ewald->enlarge( Vec3d{0.0,0.0,20.0} ); // add 20A of Vaccuum to simulated slabe (not periodic in z direction)
+        ewald->potential_of_atoms( grid.n.z, VCoul, natoms_, apos_, qs.data() );
+    }
+
+#else
+
+    void makeCoulombEwald(  int natoms_, Vec3d * apos_, Quat4d * REQs_, double* VCoul=0 ){
+        printf( "GridFF::makeCoulombEwald() ERROR: FFTW is not compiled in this version of GridFF\n" );
+        exit(0);
+    }
+
+#endif 
+
+
+
+
     void tryLoadGridFF_potentials( int natoms_, Vec3d * apos_, Quat4d * REQs_, int ntot, double*& VPaul, double*& VLond, double*& VCoul, bool bSaveNPY=false, bool bSaveXSF=false, bool bForceRecalc=false ){
         printf( "GridFF::tryLoadGridFF_potentials() natoms_=%i bUseEwald=%i bForceRecalc=%i bSaveNPY=%i bSaveXSF=%i \n", natoms_, bUseEwald, bForceRecalc, bSaveNPY, bSaveXSF );
         //bool bPBC  = false;
@@ -1026,12 +1050,13 @@ inline void addForce( const Vec3d& pos, const Quat4f& PLQ, Quat4f& fe ) const {
             _allocIfNull( VCoul, ntot );
             if(bUseEwald){
                 evalBsplineRef( natoms_, apos_, REQs_, VPaul, VLond, 0 );
-                std::vector<double> qs(natoms_); for( int i=0; i<natoms_; i++ ){ qs[i] = REQs_[i].z; }
-                //std::vector<Vec3d>  ps(natoms_); for( int i=0; i<natoms_; i++ ){  Vec3d p=apos_[i];  p.z-=grid.dCell.c.z; ps[i]=p; } // NOTE: For some unknown reason we have to shift atoms by 1 pixel in z-direction (perhaps something about grid alignement)
-                ewald->copy( grid );
-                ewald->enlarge( Vec3d{0.0,0.0,20.0} ); // add 20A of Vaccuum to simulated slabe (not periodic in z direction)
-                //ewald->potential_of_atoms( grid.n.z, VCoul, natoms_, ps.data(), qs.data() );
-                ewald->potential_of_atoms( grid.n.z, VCoul, natoms_, apos_, qs.data() );
+                // std::vector<double> qs(natoms_); for( int i=0; i<natoms_; i++ ){ qs[i] = REQs_[i].z; }
+                // //std::vector<Vec3d>  ps(natoms_); for( int i=0; i<natoms_; i++ ){  Vec3d p=apos_[i];  p.z-=grid.dCell.c.z; ps[i]=p; } // NOTE: For some unknown reason we have to shift atoms by 1 pixel in z-direction (perhaps something about grid alignement)
+                // ewald->copy( grid );
+                // ewald->enlarge( Vec3d{0.0,0.0,20.0} ); // add 20A of Vaccuum to simulated slabe (not periodic in z direction)
+                // //ewald->potential_of_atoms( grid.n.z, VCoul, natoms_, ps.data(), qs.data() );
+                // ewald->potential_of_atoms( grid.n.z, VCoul, natoms_, apos_, qs.data() );
+                makeCoulombEwald( natoms_, apos_, REQs_, VCoul );
             }else{
                 evalBsplineRef( natoms_, apos_, REQs_, VPaul, VLond, VCoul );
             }

cpp/common_SDL/SDL2OGL/GUI.h

@@ -10,6 +10,7 @@
 #include "Draw3D.h"
 
 #include "Table.h"
+#include "Interfaces.h"
 
 #include <string>
 #include <functional>
@@ -527,6 +528,59 @@ class TableView : public GUIAbstractPanel { public:
 
 };
 
+// ==============================
+//    BoundGUI & GUIPanelWatcher
+// ==============================
+
+/**
+ * @class GUIPanelWatcher
+ * @brief A class that represents a watcher for GUIPanel objects.
+ * 
+ * The GUIPanelWatcher class is used to bind a GUIPanel object with a slave object and synchronize their values.
+ * It provides methods to bind, load, and apply the values between the master and slave objects.
+ * The watcher can handle both integer and double values based on the provided flag.
+ */
+class GUIPanelWatcher{ public:
+    GUIPanel* master;  // GUIPanel object which visualize state of the slave object and can change slave's state
+    void*     slave;   // Pointer to data value to be watched (i.e. visualized and controlled by the master GUIPanel object)
+    bool bInt=false;
+    void bind    ( GUIPanel* master_, void* slave_, bool bInt_ ){ master=master_; slave=slave_; bInt=bInt_; };
+    void bindLoad( GUIPanel* master_, void* slave_, bool bInt_ ){ bind(master_,slave_,bInt_); load(); master_->redraw=true; };
+    void apply(){ if(bInt){ *(int*)slave  =  (int )master->value;   }else{ *(double*)slave = master->value;   };                     };   // change slave value  
+    void load (){ if(bInt){ master->value = *(int*)slave;           }else{ master->value   = *(double*)slave; }; master->val2text(); };   // read and visualize slave value 
+    bool check(){ if(master&&slave) if( master->redraw ){ apply(); return true; }; return false; }    // update if master should be redrawn
+};
+
+/**
+ * @class BoundGUI
+ * @brief A class representing a bound graphical user interface to some set of parameters
+ * 
+ * This class inherits from MultiPanel and BindLoader.
+ * It provides functionality for managing a GUI with multiple panels and binding data to the GUI elements.
+ */
+class BoundGUI:public MultiPanel,public BindLoader{ public:
+    bool binded=false;
+    GUIPanelWatcher* drivers=0;  // list of GUIPanelWatcher which control individual binded values, each correspond to one sub-panel of MultiPanel
+    BoundGUI(const std::string& caption, int xmin, int ymin, int xmax, int dy,int nsub):MultiPanel(caption,xmin,ymin,xmax,dy,nsub){
+        opened=false;
+    }
+    virtual void view()override{ if(opened)MultiPanel::view(); };
+    //virtual int bindLoad(void* o)=0;
+    void unbind(){
+        binded=false;
+        close();
+        //opened=false;
+        //redraw=true; tryRender();
+    }
+    bool check(){
+        if(!binded) return false;
+        bool bChanged=false;
+        for(int i=0; i<nsubs; i++){
+            bChanged |= drivers[i].check();
+        }
+        return bChanged;
+    }
+};
 
 // ==============================
 //    class GUI

doc/Chebyshev_Acceleration.md

@@ -0,0 +1,190 @@
+# Chebyshev-Accelerated Jacobi Solver for Projective Dynamics
+
+## 1. Introduction
+
+This document provides a complete implementation guide for the Chebyshev-accelerated Jacobi solver as described in [Wang 2015] (Section 2 of the source document). The method combines Jacobi iterations with Chebyshev semi-iterative acceleration for fast convergence in projective dynamics simulations.
+
+## 2. Core Algorithm
+
+## 2.1 Chebyshev Acceleration Formula
+
+From Section 2.1, the definitive Chebyshev acceleration formula is:
+
+
+
+lets express this in terms of displacement vectors $d_k$ by which the positions $q_k$ are changed in each iteration:
+
+$ q_{k+1} = ω_{k+1} ( γ {\hat d}_{k+1} + d_k ) + q_{k−1} $
+
+$\hat d_{k+1} = q̂_{k+1}- q_k$
+$d_{k+1} = q_{k+1}- q_k$
+$d_k     = q_k- q_{k-1}$
+
+
+Elimineate $q_{k-1}$ by using 
+$ q_{k-1} = q_k - d_k $
+
+$ q_{k+1} = ω ( γ ( q̂_{k+1} - q_k ) + q_k - q_k - d_k ) + q_k - d_k $
+
+
+$ q_{k+1} = ω  γ  q̂_{k+1} - ω  γ q_k  + ω q_k - ω q_k - ω d_k  + q_k - d_k $
+
+$ q_{k+1} = ω  γ  q̂_{k+1} + ( 1 - ω  γ ) q_k         - (1+ω) d_k $
+
+
+
+
+
+$ q_{k+1} = ω_{k+1} ( γ ( q̂_{k+1} - q_k ) + q_k - q_{k−1} ) + q_{k−1} $
+$ q_{k+1} = ω_{k+1}  γ ( q̂_{k+1} - q_k ) + ω_{k+1} q_k - ω_{k+1} q_{k−1} + q_{k−1} $
+$ q_{k+1} = ω_{k+1}    q̂_{k+1} -  ω_{k+1} ( 1  - γ ) q_k + (1 - ω_{k+1}) q_{k−1} $
+
+$ q_{k+1} = ω ( γ ( q̂_{k+1} - q_k ) + q_k - q_{k−1} ) + q_{k−1} $
+$ q_{k+1} = ω   γ ( q̂_{k+1} - q_k ) + ω q_k - ω q_{k−1} + q_{k−1} $
+$ q_{k+1} = ω   q̂_{k+1}     -  ω ( 1  - γ ) q_k + (1 - ω) q_{k−1} $
+
+where:
+ - $x_{new}$​ is the result of one Jacobi iteration
+ - $x_{k}$​ is the previous position
+ - $ω_{k+1}$​ is the Chebyshev weight
+
+### 2.2 Chebyshev Weights Calculation
+
+The weights are calculated as:
+
+```python
+def calculate_omega(k, rho, prev_omega):
+    if k < DELAY_START:  # typically 10 iterations
+        return 1.0
+    elif k == DELAY_START:
+        return 2.0 / (2.0 - rho**2)
+    else:
+        return 4.0 / (4.0 - rho**2 * prev_omega)
+```
+
+To estimate the spectral radius (ρ), we use power iteration:
+
+
+```python
+def estimate_spectral_radius(self, A, x0):
+        """Estimate spectral radius using power iteration"""
+        x = x0.copy()
+        x_prev = x0.copy()
+
+        # Apply matrix twice
+        x = A @ x
+        x_next = A @ x
+
+        # Estimate using Rayleigh quotient
+        rho = np.sqrt(abs(np.dot(x_next, x) / np.dot(x, x)))
+        return min(rho, 0.9992)  # Cap at 0.9992 as suggested in Section 2.2.2
+```
+
+
+### 2.3 Complete Implementation
+
+```python
+
+def solve(self, A, b, x0, max_iter=400, tol=1e-6):
+    """
+    Solve Ax = b using Chebyshev-accelerated Jacobi
+    Args:
+        A: System matrix
+        b: Right-hand side vector
+        x0: Initial guess
+    """
+    # Initialize
+    x_prev = x0.copy()
+    x_curr = x0.copy()
+
+    # Estimate spectral radius
+    rho = estimate_spectral_radius(A, x0)
+
+    # Initialize omega
+    omega = 1.0
+
+    for k in range(max_iter):
+        # Regular Jacobi iteration
+        D_inv = 1.0 / np.diag(A)
+        x_new = D_inv * (b - (A @ x_curr - np.diag(A) * x_curr))
+
+        # Calculate new Chebyshev weight
+        omega_next = self.calculate_omega(k, rho, omega)
+
+        # Chebyshev acceleration
+        x_accel = omega_next * (x_new - x_prev) + x_prev
+
+        # Update positions
+        x_prev = x_curr
+        x_curr = x_accel
+        omega = omega_next
+
+        # Check convergence
+        residual = np.linalg.norm(b - A @ x_curr)
+        if residual < tol:
+            break
+
+    return x_curr
+```
+
+## 3. Important Implementation Notes
+
+### 3.1 Spectral Radius (ρ)
+
+From Section 2.2.2:
+
+ - Use ρ ≈ 0.9992-0.9999 for most cases
+ - Underestimating ρ is better than overestimating
+ - Can be estimated using power iteration as shown above
+
+### 3.2 Delayed Start
+
+From Section 2.2.2 of the paper:
+
+ - Delay Chebyshev acceleration for first S iterations (typically 10)
+ - Set ω₁ = ω₂ = ... = ωₛ = 1 for these iterations
+ - Helps prevent oscillation/divergence in early iterations
+
+### 3.3 Convergence Properties
+
+From Section 2.2.2:
+
+ - Method converges faster in first few iterations
+ - Convergence rate drops after initial phase
+ - Oscillation indicates ρ is too large
+ - Divergence is rare if ρ is properly tuned
+
+4. Usage Example
+
+```python
+# Create solver
+solver = ChebyshevAcceleratedSolver(delay_start=10)
+
+# Setup system
+A = create_system_matrix()  # Your system matrix
+b = create_rhs_vector()     # Your right-hand side
+x0 = np.zeros(n)           # Initial guess
+
+# Solve system
+x = solver.solve(A, b, x0, max_iter=400, tol=1e-6)
+```
+
+## 5. Advantages and Limitations
+
+### 5.1 Advantages (Section 2.4.2)
+
+ - Highly parallelizable
+ - Small memory footprint
+ - No need for expensive linear algebra libraries
+ - More efficient than CG on GPU
+ - Robust to small system changes
+
+### 5.2 Limitations
+
+ - May require underrelaxation for stability in some cases
+ - Convergence affected by mesh quality
+ - May need smaller timesteps for highly stiff systems
+
+## References
+
+Wang, H. (2015). A Chebyshev Semi-iterative Approach for Accelerating Projective and Position-based Dynamics. ACM Trans. Graph. (SIGGRAPH Asia) 34, 6, Article 246.