Discussion :

[Kafkana1]

Bonjour,

j'ai un problème de compilation avec des spécialisations template. Je bosse sur un code de dynamique moléculaire et j'ai des espaces en 2 et 3 dimensions. Le probleme, c'est que je n'ai pas les même type pour ma variable de rotation (double en 2D et tinyvector en 3D) et quand je veux compiler mon code en 2D, le compilateur essaie aussi de compiler la version 3D de mon .cpp.

erreur: no matching function for call to ‘mols::PatchyParticleHybrid<double, 2>::rotate(mols::PatchyParticleHybrid<double, 2>::Position)’

J'aimerais donc savoir si il y a un moyen de contourner le probleme.

.cpp

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
template<int D>
void addMonomer(Config& config){}
 
template<>
void addMonomer<2>(Config& config) 
{
     //autres trucs 
p.rotate(typename Config::AtomType::positionType (M_PI)); 
 
}
 
template<>
void addMonomer<3>(Config& config) 
{
     //autres trucs 
p.rotate(typename Config::AtomType::Position (0,0,M_PI)); //<--- probleme ici quand D = 2
 
}

.hpp

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
 template <typename pType, int D>
      class Particle : public Atom<pType, D>
      {};
 
      template <typename pType>
      class Particle<pType, 2> : public Atom<pType, 2>
      {
	    typedef pType positionType;
 
	    inline void rotate(positionType angle)
	    {
		//des trucs
	    }
      }
 
 
      template <typename pType>
      class Particle<pType, 3> : public Atom<pType, 3>
      {
	    typedef blitz::TinyVector<positionType, 3> Position;
 
	    inline void rotate(Position angle)
	    {
		//des trucs
	    }
      }

merci à vous ^^

[darkman19320]

Salut,

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

typedef blitz::TinyVector<positionType, 3> Position;

Où est-ce que positionType est défini dans ta spécialisation de classe:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
template <typename pType>
class Particle<pType, 3> : public Atom<pType, 3>

En plus, ta méthode rotate est privée donc difficilement utilisable.

Deuxième point, pourquoi faire un typedef dans une spécialisation de template alors que tu pourrais avoir un truc du style:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
 
template <typename pType, int D>
class Particle : public Atom<pType, D>
{
public:
    void rotate(const pType& angle)
    {
       //des trucs du genre si D = 2 appelle d'une méthode 2D, sinon 3D.
    }
};
 
typedef Particule<double,2> Particule2D;
typedef Particule<blitz::TinyVector<double, 3> > Particule3D;

Dernière chose, les templates doivent être définie dans le fichier d'entête.

Bon courage!

[Kafkana1]

Merci pour ta réponse ^^

En fait, je bosse sur un gros programme utilisé et modifié par plusieurs personnes. Seulement, on a un "coeur" commun. Je me suis donc adapté au code et les typedef servent dans d'autres parties du code, auxquelles je ne peux pas toucher.

j'ai été un peu vite la derniere fois, je m'en excuse, je vais mettre les parties de code completes

.hpp

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
      template <typename pType, int D>
      class Particle : public Atom<pType, D>
      {};
 
      /* Class dédiée aux particles 2D. C'est une sous classe d'atome dans laquelle on rajoute :
 
	 - une vitesse angulaire
	 - un moment d'intertie
	 - un angle de rotation
	 - un couple
 
      */
 
      template <typename pType>
      class Particle<pType, 2> : public Atom<pType, 2>
      {
 
	    typedef Atom<pType, 2> Base;
 
      public:
	    typedef pType Scalar;
	    typedef pType positionType;
	    typedef pType velocityType;
            typedef blitz::TinyVector<positionType, 2> Position;
            typedef blitz::TinyVector<velocityType, 2> Velocity;
            typedef blitz::TinyVector<velocityType, 2> Force;
 
	    Particle() : Atom<pType, 2>(), angular_velocity_(0.0), rotational_angle_(0.0),  moment_of_inertia_(1.0), torque_(0.0)
	    {}
 
            Particle(Scalar r, Scalar density = 1.0) : Atom<pType, 2>(r) {
                  Base::mass()              = density * M_PI * Base::radius() * Base::radius();
                  moment_of_inertia() = 0.5 * Base::mass() * Base::radius() * Base::radius();		  
            }
 
            Particle(Scalar r, twister& tw, double T, Scalar density = 1.0) : Atom<pType, 2>(r) {
                  Base::mass()              = density * M_PI * Base::radius() * Base::radius();
                  moment_of_inertia() = 0.5 * Base::mass() * Base::radius() * Base::radius();
                  rotational_angle_         = tw.drand() * 2 * M_PI;
 
                  angular_velocity()  = tw.grand()*sqrt(T/moment_of_inertia());
                  for(int d = 0; d<2; d++)
                        Base::velocity()[d] = tw.grand()*sqrt(T/Base::mass());
            }
 
	    inline       velocityType& angular_velocity()            { return angular_velocity_; }
	    inline const velocityType& angular_velocity()      const { return angular_velocity_; }
	    inline       Scalar&       moment_of_inertia()           { return moment_of_inertia_;}
	    inline const Scalar&       moment_of_inertia()     const { return moment_of_inertia_;}
	    inline const positionType& rotational_angle()      const { return rotational_angle_; }
	    inline       Scalar&       torque()                      { return torque_;           }
	    inline const Scalar&       torque()                const { return torque_;           }
 
 
	    inline void clearForces() {  Base::clearForces(); 
		  torque_ = 0.0;      }
 
 
	    template <class PairType>
	    inline void collect1(const PairType& pair) { Base::force_  += pair.fij();
		  torque_ += pair.Ci();} 
 
 
 
	    template <class PairType>
	    inline void collect2(const PairType& pair) { Base::force_  -= pair.fij(); 
		  torque_ -= pair.Cj(); }
 
 
	    // on a décidé de ne pas modifier les rotations hors de la classe	    
	    inline void rotate(Scalar angle)
	    {
		  rotational_angle_ = rotational_angle_ + angle;  
	    }
 
 
	    void readProperties(ibstream& in)
	    {
		  Base::readProperties(in);
		  if(in.mode()&iomode::output_masses) in.read(moment_of_inertia_);
	    }
 
	    void writeProperties(obstream& out) const
	    {
		  Base::writeProperties(out);
		  if(out.mode()&iomode::output_masses) out.write(moment_of_inertia_);
	    }
 
 
	    void readData(ibstream& in) 
	    {
		  Base::readData(in);
		  in.read(rotational_angle_);
		  if(in.mode()&iomode::output_velocities) in.read(angular_velocity_);
	    }
 
	    void writeData(obstream& out) const
	    {
		  Base::writeData(out);
		  out.write(rotational_angle_);
		  if(out.mode()&iomode::output_velocities) out.write(angular_velocity_);
	    }
 
 
 
      protected: 
 
	    velocityType          angular_velocity_;
	    positionType          rotational_angle_;
	    Scalar                moment_of_inertia_;
	    Scalar                torque_;
      };
 
 
 
 
      /* Class dédiée aux particles 3D. C'est une sous classe d'atome dans laquelle on rajoute :
 
	 - une vitesse angulaire
	 - un moment d'intertie
	 - un angle de rotation
	 - un couple
 
      */
 
      template <typename pType>
      class Particle<pType, 3> : public Atom<pType, 3>
      {
 
	    typedef Atom<pType, 3> Base;
 
      public:
	    typedef pType Scalar;
	    typedef pType positionType;
	    typedef pType velocityType;
            typedef blitz::TinyVector<positionType, 3> Position;
            typedef blitz::TinyVector<velocityType, 3> Velocity;
            typedef blitz::TinyVector<velocityType, 3> Force;
 
 
	    Particle() : Atom<pType, 3>(), angular_velocity_(0.0),  moment_of_inertia_(1.0), torque_(0.0)
	    {}
 
	    Particle(Scalar r, Scalar density = 1.0) : Atom<pType, 3>(r), angular_velocity_(0.0), torque_(0.0) {
                  Base::mass()        = ( 4.0 / 3.0 ) * density * M_PI * Base::radius() * Base::radius() * Base::radius();
                  moment_of_inertia() = 0.4 * Base::mass() * Base::radius() * Base::radius();
		  q() = 1.0, 0.0, 0.0, 0.0;		  
            }
 
            Particle(Scalar r, twister& tw, double T, Scalar density = 1.0) : Atom<pType, 3>(r) {
                  Base::mass()              = ( 4.0 / 3.0 ) * density * M_PI * Base::radius() * Base::radius() * Base::radius();
                  moment_of_inertia()       = 0.4 * Base::mass() * Base::radius() * Base::radius();
		  double phi   = tw.drand() * 2 * M_PI;
		  double theta = tw.drand() * 2 * M_PI;
		  double psi   = tw.drand() * 2 * M_PI;
		  q() = cos(0.5 * theta) * cos(0.5 * (phi+psi)), sin(0.5 * theta) * cos(0.5 * (phi-psi)), sin(0.5 * theta) * sin(0.5 * (phi-psi)), cos(0.5 * theta) * sin(0.5 * (phi+psi));
 
                  for(int d = 0; d<3; d++) {
                        angular_velocity()[d] = tw.grand()*sqrt(T/moment_of_inertia());
                        Base::velocity()[d]   = tw.grand()*sqrt(T/Base::mass());
                  }
            }
 
	    inline       Velocity&     angular_velocity()                   { return angular_velocity_;  }
	    inline const Velocity&     angular_velocity()             const { return angular_velocity_;  }
	    inline       Scalar&       moment_of_inertia()                  { return moment_of_inertia_; }
	    inline const Scalar&       moment_of_inertia()            const { return moment_of_inertia_; }
	    inline       Force&        torque()                             { return torque_;            }
	    inline const Force&        torque()                       const { return torque_;            }
	    inline const Scalar*                              A()     const { return &A_[0][0];          }
	    inline       Scalar*                              A()           { return &A_[0][0];          }
	    inline const blitz::TinyVector<Scalar, 4>&        q()     const { return q_;                 }
	    inline       blitz::TinyVector<Scalar, 4>&        q()           { return q_;                 }
 
 
	    inline void clearForces() {  Base::clearForces(); 
		  torque_ = 0.0;      }
 
 
	    template <class PairType>
	    inline void collect1(const PairType& pair) { Base::force_  += pair.fij();
		  //std::cout << "test fij struc" << Base::force_  << std::endl;
		  torque_ += pair.Ci();} 
 
 
 
	    template <class PairType>
	    inline void collect2(const PairType& pair) { Base::force_  -= pair.fij(); 
		  torque_ -= pair.Cj(); }
 
 
	    inline void rotate(Position angle)
	    {
		  angle *= 0.5;
		  q_[0] -=  q_[1] * angle[0] + q_[2] * angle[1] + q_[3] * angle[2];
		  q_[1] +=  q_[0] * angle[0] - q_[3] * angle[1] + q_[2] * angle[2];
		  q_[2] +=  q_[3] * angle[0] + q_[0] * angle[1] - q_[1] * angle[2];
		  q_[3] += -q_[2] * angle[0] + q_[1] * angle[1] + q_[0] * angle[2];
		  q_ /= sqrt(blitz::sum(q_*q_));
	    }					     
 
 
	    void readProperties(ibstream& in)
	    {
		  Base::readProperties(in);
		  if(in.mode()&iomode::output_masses) in.read(moment_of_inertia_);
	    }
 
	    void writeProperties(obstream& out) const
	    {
		  Base::writeProperties(out);
		  if(out.mode()&iomode::output_masses) out.write(moment_of_inertia_);
	    }
 
 
	    void readData(ibstream& in) 
	    {
		  Base::readData(in);
		  pType phi, theta, psi;
		  in.read(phi);
		  in.read(theta);
		  in.read(psi);
		  q() = cos(0.5 * theta) * cos(0.5 * (phi+psi)), sin(0.5 * theta) * cos(0.5 * (phi-psi)), sin(0.5 * theta) * sin(0.5 * (phi-psi)), cos(0.5 * theta) * sin(0.5 * (phi+psi));
 
		  if(in.mode()&iomode::output_velocities) in.read(angular_velocity_);
	    }
 
	    void writeData(obstream& out) const
	    {		 
		  Base::writeData(out);
		  pType q00 = q_[0]*q_[0];
		  pType q11 = q_[1]*q_[1];
		  pType q22 = q_[2]*q_[2];
		  pType q33 = q_[3]*q_[3];
		  pType q01 = q_[1]*q_[0];
		  pType q02 = q_[0]*q_[2];
		  pType q03 = q_[0]*q_[3];
		  pType q12 = q_[1]*q_[2];
		  pType q13 = q_[1]*q_[3];
		  pType q23 = q_[2]*q_[3];
		  pType phi  = atan2(2*(q13+q02), -2*(q23-q01));
		  pType theta= acos((q00-q11-q22+q33));
		  pType psi  = atan2(2*(q13-q02), 2*(q23+q01));
		  out.write(phi);
		  out.write(theta);
		  out.write(psi);
		  if(out.mode()&iomode::output_velocities) out.write(angular_velocity_);
	    }
 
	    void updateA() 
	    {
		  pType q00 = q_[0]*q_[0];
		  pType q11 = q_[1]*q_[1];
		  pType q22 = q_[2]*q_[2];
		  pType q33 = q_[3]*q_[3];
		  pType q01 = q_[1]*q_[0];
		  pType q02 = q_[0]*q_[2];
		  pType q03 = q_[0]*q_[3];
		  pType q12 = q_[1]*q_[2];
		  pType q13 = q_[1]*q_[3];
		  pType q23 = q_[2]*q_[3];
 
		  A_[0][0] =    q00 + q11 - q22 - q33;
		  A_[0][1] =    2 * (q12 + q03);
		  A_[0][2] =    2 * (q13 - q02);
 
		  A_[1][0] =    2 * (q12 - q03);
		  A_[1][1] =    q00 - q11 + q22 - q33;
		  A_[1][2] =    2 * (q01+q23);
 
		  A_[2][0] =   2 * (q13 + q02);
		  A_[2][1] =   2 * (- q01 + q23);
		  A_[2][2] =   q00 - q11 - q22 + q33;
	    }
 
      protected: 
 
 
	    Velocity                                               angular_velocity_;
	    Scalar                                                 moment_of_inertia_;
	    Force                                                  torque_;
	    pType                                                  A_[4][4];  
	    blitz::TinyVector<Scalar, 4>                           q_;
      };
 
 
      template <typename pType, int D> struct recursion_level< Particle<pType,D> > { typedef atom_tag type; };

et mon .cpp

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
twister tw;
double T=0.05, eta = 0.01, dt=0.001, gdot =0.0, radOne=0.495, radTwo=0.5, pf;
int N, M;
 
template<const int D>
void addMonomer(Config& config){}
 
template<>
void addMonomer<1>(Config& config) 
{
      std::cerr << "Dimension = " << dim << " not implemented " << std::endl;
}
 
template<>
void addMonomer<2>(Config& config) 
{
      int n = config.atoms().size();
      double e = energy(config);
      if((n+M)%M == 0) 
	    {
		  PatchyParticleHybrid<Scalar, dim> p(1, radOne, tw, T);
		  p.position() = tw.drandVec<Scalar, dim>()*config.cell().side();
		  config.atoms().push_back(p);
	    }
 
      else 
	    {
		  PatchyParticleHybrid<Scalar, dim> p = config.atoms()[n-1];
		  if((n+1)%M != 0) 
			{
			      p.radius() = radTwo;
			      p.mass() =  M_PI *p.radius()*p.radius();
			      p.moment_of_inertia() = 0.5*p.mass()*p.radius()*p.radius();
			      p.setPatches(2);
			      p.position() += 2.08*p.patch(0);
			}
		  else 
			{
			      p.radius() = radOne;
			      p.mass() =  M_PI *p.radius()*p.radius();
			      p.moment_of_inertia() = 0.5*p.mass()*p.radius()*p.radius();
			      p.setPatches(1);
			      p.position() += 2.08*p.patch(0);
			      p.rotate(typename Config::AtomType::positionType (M_PI));
			}
 
 
		  config.atoms().push_back(p);
	    }
 
      config.reset();
      //cout << n << " " << energy(config)-e << endl;
      if(energy(config) > e) 
	    {
		  config.atoms().resize(n); // nothing done
		  config.reset();
	    }
}
 
template<>
void addMonomer<3>(Config& config) 
{
 
      int n = config.atoms().size();
      double e = energy(config);
      if((n+M)%M == 0) 
	    {
		  PatchyParticleHybrid<Scalar, dim> p(1, radOne, tw, T);
		  p.position() = tw.drandVec<Scalar, dim>()*config.cell().side();
		  config.atoms().push_back(p);
	    }
 
      else 
	    {
		  PatchyParticleHybrid<Scalar, dim> p = config.atoms()[n-1];
		  if((n+1)%M != 0) 
			{
			      p.radius() = radTwo;
			      p.mass() =  ( 4.0 / 3.0 ) * M_PI *p.radius()*p.radius()*p.radius();
			      p.moment_of_inertia() = 0.4*p.mass()*p.radius()*p.radius();
			      p.setPatches(2);
			      p.position() += 2.08*p.patch(0);
			}
		  else 
			{
			      p.radius() = radOne;
			      p.mass() =  ( 4.0 / 3.0 ) * M_PI *p.radius()*p.radius()*p.radius();
			      p.moment_of_inertia() = 0.4*p.mass()*p.radius()*p.radius();
			      p.setPatches(1);
			      p.position() += 2.08*p.patch(0);
			      p.rotate(typename Config::AtomType::Position (0,0,M_PI)); 
			}	
 
		  config.atoms().push_back(p);
	    }
 
      config.reset();
      //cout << n << " " << energy(config)-e << endl;
      if(energy(config) > e) 
	    {
		  config.atoms().resize(n); // nothing done
		  config.reset();
	    }
}
 
 
int main(int argc, char** argv) {
//pas mal de choses
addMonomer<dim>(config);
//pas mal d'autres choses

il est vrai que j'aurais du mettre ma fonction addmonomer dans un hpp à part.

Le problème, c'est que si je supprime entierement la fonction addmonomer<3> (la dimension étant choisi a la compilation avec l'option -DSPACEDIM) mon code marche très bien en 2D. (il marche de toute façon en 3D car blitz::tinyvector accepte les doubles en entrée). J'aimerais donc savoir si il existe un moyen d'empêcher le compilateur d'essayer d'acceder à la ligne

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

p.rotate(typename Config::AtomType::Position (0,0,M_PI));

qui ne me sert qu'en compilation 3D, quand je compile en 2D?

[darkman19320]

Vu que tu utilises une condition préprocesseur, tu vas pouvoir activer/désactiver des parties de ton code.
Avec la définition SPACEDIM, tu as tous les éléments pour le faire.

Voici un exemple:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
#if SPACEDIM == 1
template<>
void addMonomer<1>(Config& config) 
{
      std::cerr << "Dimension = " << dim << " not implemented " << std::endl;
}
#elif SPACEDIM == 2
template<>
void addMonomer<2>(Config& config) 
{
//...
}
#elif SPACEDIM == 3
template<>
void addMonomer<3>(Config& config) 
{
//...
}
#endif

Mais dis-toi bien que ça va rendre ton code pas super lisible si tu as de plus en plus de fonctions de ce genre.

[Kafkana1]

super ça marche!

merci pour l'aide et les conseils !