FhSim  3.1.0
Marine systems simulation
Loading...
Searching...
No Matches
PotentialData.h
1#ifndef _CPOTENTIAL_DATA_H_
2#define _CPOTENTIAL_DATA_H_
3
4#include <iostream>
5#include <iomanip>
6#include <fstream>
7#include <algorithm>
8#include <string>
9#include <vector>
10#include <complex>
11#include <limits>
12#include <map>
13
14#include <Eigen/Eigen>
15#include <Eigen/Eigenvalues>
16#include <Eigen/Core>
17
18
19#include "../HullData/HullData.h"
20
21#include "./Wave/FroudeKrylov.h"
22#include "./Wave/WaveDrift.h"
23
24#include "Identification/FreqDomainIdent.h"
25
26
27namespace Ship{
28 class PotentialData
29 {
30 public:
31 EIGEN_MAKE_ALIGNED_OPERATOR_NEW
32
33 struct TimeDomainRadiation
34 {
35 EIGEN_MAKE_ALIGNED_OPERATOR_NEW
36
37 TimeDomainRadiation()
38 {
39 Speed = -10;
40 Ainf = Eigen::Matrix< double, 6, 6 >::Zero( );
41
42 for( int i=0; i <6; i++)
43 for( int j=0; j <6; j++)
44 {
45 A[i][j] = Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic >::Zero( 1, 1 );;
46 B[i][j] = Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic >::Zero( 1, 1 );;
47 C[i][j] = Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic >::Zero( 1, 1 );;
48 D[i][j] = Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic >::Zero( 1, 1 );;
49 }
50 }
51
52 // Ainf
53 Eigen::Matrix<double, 6, 6> Ainf;
54
55 double Speed;
56
57 // companion system matrices
58 Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> A[6][6];
59 Eigen::Matrix<double, Eigen::Dynamic, 1 > B[6][6];
60 Eigen::Matrix<double, 1 , Eigen::Dynamic> C[6][6];
61 Eigen::Matrix<double, Eigen::Dynamic, 1 > D[6][6];
62 double Direction;
63 };
64
65 virtual ~PotentialData()
66 {
67 delete Speed;
68 delete Heading;
69 delete Frequency;
70 }
71
72 virtual void Read( const std::string &fname, int flipPhaseForOSC =0) = 0;
73 virtual Eigen::Matrix<double,6,6> GetHToO() = 0;
74
76 const std::vector<std::string>& Warnings() const { return m_warnings; }
77 void ClearWarnings() { m_warnings.clear(); }
78
79 virtual void SetOriginWorigin() = 0;
80 virtual void SetOriginGorigin() = 0;
81 virtual void SetOrigin( Eigen::Matrix<double,3,1> &r ) = 0;
82
83 Eigen::Matrix<double,6,6> GetMassRigidBody()
84 {
85 return GetHToO().transpose()*Mrbeig*GetHToO();
86 }
87
88 Eigen::Matrix<double,6,6> GetAinf(unsigned int speedno=0)
89 {
90 return RadiationApproximation[speedno].Ainf;
91 }
92
93 Eigen::Matrix<double,6,6> GetMassTotal(unsigned int speedno=0)
94 {
95 return RadiationApproximation[speedno].Ainf+GetHToO().transpose()*Mrbeig*GetHToO();
96 }
97
98 Eigen::Matrix<double,6,6> GetAw(unsigned int FreqNo, unsigned int SpeedNo=0)
99 {
100 return GetHToO().transpose()*A[SpeedNo][FreqNo]*GetHToO();
101 }
102
103 Eigen::Matrix<double,6,6> GetBw(unsigned int FreqNo, unsigned int SpeedNo=0)
104 {
105 return GetHToO().transpose()*B[SpeedNo][FreqNo]*GetHToO();
106 }
107
108 virtual Eigen::Matrix<double,6,6> GetCw(unsigned int FreqNo,unsigned int SpeedNo=0)
109 {
110 return GetHToO().transpose()*C[SpeedNo][FreqNo]*GetHToO();
111 }
112
113 Eigen::Matrix<double,1,Eigen::Dynamic> GetAij(unsigned int DOF1, unsigned int DOF2, unsigned int SpeedNo=0)
114 {
115 Eigen::Matrix<double,1,Eigen::Dynamic> A = Eigen::Matrix<double,1,Eigen::Dynamic>::Zero(1,NoFrequency);
116
117 for( unsigned int i=0; i < NoFrequency; i++)
118 {
119 Eigen::Matrix<double,6,6> Aw = GetAw(i,SpeedNo);
120 A(i) = Aw(DOF1,DOF2);
121 }
122 return A;
123 }
124
125 Eigen::Matrix<double,1,Eigen::Dynamic> GetBij(unsigned int DOF1, unsigned int DOF2, unsigned int SpeedNo=0)
126 {
127 Eigen::Matrix<double,1,Eigen::Dynamic> B = Eigen::Matrix<double,1,Eigen::Dynamic>::Zero(1,NoFrequency);
128
129 for( unsigned int i=0; i < NoFrequency; i++)
130 {
131 Eigen::Matrix<double,6,6> Bw = GetBw(i,SpeedNo);
132 B(i) = Bw(DOF1,DOF2);
133 }
134 return B;
135 }
136
137 Eigen::Matrix<double,1,Eigen::Dynamic> GetCij(unsigned int DOF1, unsigned int DOF2, unsigned int SpeedNo=0)
138 {
139 Eigen::Matrix<double,1,Eigen::Dynamic> C = Eigen::Matrix<double,1,Eigen::Dynamic>::Zero(1,NoFrequency);
140
141 for( unsigned int i=0; i < NoFrequency; i++)
142 {
143 Eigen::Matrix<double,6,6> Cw = GetCw(i,SpeedNo);
144 C(i) = Cw(DOF1,DOF2);
145 }
146 return C;
147 }
148
149 Eigen::Matrix<double,1,Eigen::Dynamic> GetSpeed()
150 {
151 Eigen::Matrix<double,1,Eigen::Dynamic> speed = Eigen::Matrix<double,1,Eigen::Dynamic>::Zero(1,NoSpeed);
152 for(unsigned int i=0; i < NoSpeed; i++ ) speed(i) = Speed[i];
153 return speed;
154 }
155
156 Eigen::Matrix<double,1,Eigen::Dynamic> GetHeading()
157 {
158 Eigen::Matrix<double,1,Eigen::Dynamic> head = Eigen::Matrix<double,1,Eigen::Dynamic>::Zero(1,NoHeading);
159 for(unsigned int i=0; i < NoHeading; i++ ) head(i) = Heading[i];
160 return head;
161 }
162
163 Eigen::Matrix<double,1,Eigen::Dynamic> GetFrequency(int SpeedNo=0)
164 {
165 Eigen::Matrix<double,1,Eigen::Dynamic> freq = Eigen::Matrix<double,1,Eigen::Dynamic>::Zero(1,NoFrequency);
166 for(unsigned int i=0; i < NoFrequency; i++ ) freq(i) = Frequency[SpeedNo][i];
167 return freq;
168 }
169
170 virtual Eigen::Matrix<double,3,1> GetRg()
171 {
172 return rg;
173 }
174
175 virtual Eigen::Matrix<double,3,1> GetRo()
176 {
177 return ro;
178 }
179
180 Eigen::Matrix<double,6,1> GetFirstOrderForce( const Eigen::Matrix<double,3,1> &Kxy, double ShipSpeed,
181 double WaveDirection, double WaveAmplitude, double WaveFrequency, double WaveNumber, double WavePhase,
182 double Time )
183 {
184 double Fptr[6];
185
186 double EncounterFrequency = WaveFrequency-WaveNumber*ShipSpeed;
187 double Direction = EncounterFrequency/std::abs(EncounterFrequency);
188 EncounterFrequency = std::abs(EncounterFrequency);
189
190 //std::cout << "T: " << Time << "\t";
191
192 for( int DOF=0; DOF < 6; DOF++)
193 {
194
195 std::complex<double> RAO = FroudeFrylov->GetRAOValue(DOF,WaveDirection,std::abs(ShipSpeed),WaveFrequency);
196
197
198 double Amplitude = std::abs(RAO);
199 double Phase = atan2( RAO.imag(), RAO.real() );
200
201 //std::cout << DOF << " " <<RAO.real() << " " << RAO.imag() << "i" << " , Lag:" << 1/WaveFrequency*(Kxy(0) - Kxy(1) - Phase - WavePhase) << "s ";
202
203 double ThisWaveHeight = WaveAmplitude*sin(EncounterFrequency *Time -Kxy(0) -Kxy(1) + Phase - WavePhase );
204 Fptr[DOF] = Direction*Amplitude*ThisWaveHeight;
205 }
206 //std::cout << std::endl;
207 Eigen::Map<Eigen::Matrix<double,6,1>,Eigen::Unaligned> tmp(Fptr);
208 return GetHToO().transpose()*tmp;
209 }
210
211 Eigen::Matrix<double,6,1> GetSecondOrderForce(Eigen::Matrix<double,3,1> ShipPosition, double ShipSpeed,
212 double WaveDirection, double WaveHeight, double WaveFrequency, double WaveNumber, double WavePhase,
213 double Time )
214 {
215 double Fptr[6] = {}; // zero without a wave drift data set (owner ruling D42)
216 for( int DOF=0; DOF < 6; DOF++)
217 {
218 if( WaveDrift )
219 {
220 std::complex<double> Amplitude = WaveDrift->GetRAOValue(DOF,ShipSpeed,WaveFrequency,WaveDirection);
221 Fptr[DOF] = std::pow(WaveHeight,2)*Amplitude.real();
222 }
223 }
224 Eigen::Map<Eigen::Matrix<double,6,1>,Eigen::Unaligned> tmp(Fptr);
225 return GetHToO().transpose()*tmp;
226 }
227
228 void GenerateZeroSpeedTimeDomainApproximation()
229 {
230 //std::cout << "Generating zero speed time domain state radiation model"<<std::endl;
231 if( !RadiationApproximation.empty() ) RadiationApproximation.clear();
232 if( ZeroSpeedIdentificationFromCache() ){
233 //std::cout << "using model from cache" << std::endl;
234 return;
235 }else{
236 //std::cout << "identifying new model" << std::endl;
237 }
238
239 FreqDomainIdent Identification;
240 Identification.SetPotentialData(this);
241 Identification.IdentificationMethod(FreqDomainIdent::VECTOR_FITTING);
242
243 TimeDomainRadiation zeroSpeedSystem;
244
245
246 for( int i=0; i < 6; i++)
247 {
248
249 for( int j=0; j <= i; j++)
250 {
251 if( (i == 1 || i == 3 || i == 5) && (j==0 || j==2 || j==4) ) continue;
252 if( (j == 1 || j == 3 || j == 5) && (i==0 || i==2 || i==4) ) continue;
253
254 std::cout << "Identifying (i/j) = (" << i << "/" << j << ") ";
255 fflush(0);
256 if( is2Ddata() && (i!=0 && j!=0)){
257 FreqDomainIdent::StateSpace SS = Identification.IdentifySS(i,j);
258 zeroSpeedSystem.Ainf(i,j) = Identification.GetAinf(i,j);
259 zeroSpeedSystem.A[i][j] = SS.A;
260 zeroSpeedSystem.B[i][j] = SS.B;
261 zeroSpeedSystem.C[i][j] = SS.C;
262 zeroSpeedSystem.D[i][j] = SS.D;
263
264 zeroSpeedSystem.Ainf(j,i) = zeroSpeedSystem.Ainf(i,j);
265 zeroSpeedSystem.A[j][i] = zeroSpeedSystem.A[i][j];
266 zeroSpeedSystem.B[j][i] = zeroSpeedSystem.B[i][j];
267 zeroSpeedSystem.C[j][i] = zeroSpeedSystem.C[i][j];
268 zeroSpeedSystem.D[j][i] = zeroSpeedSystem.D[i][j];
269 }else{
270 std::cout << " Zero data" << std::endl;
271 zeroSpeedSystem.Ainf(i,j) = GetAij(i,j,0)[0];
272 zeroSpeedSystem.A[i][j] = Eigen::Matrix<double,1,1>::Zero();
273 zeroSpeedSystem.B[i][j] = Eigen::Matrix<double,1,1>::Zero();
274 zeroSpeedSystem.C[i][j] = Eigen::Matrix<double,1,1>::Zero();
275 zeroSpeedSystem.D[i][j] = Eigen::Matrix<double,1,1>::Zero();
276
277 zeroSpeedSystem.Ainf(j,i) = GetAij(i,j,0)[0];
278 zeroSpeedSystem.A[j][i] = Eigen::Matrix<double,1,1>::Zero();
279 zeroSpeedSystem.B[j][i] = Eigen::Matrix<double,1,1>::Zero();
280 zeroSpeedSystem.C[j][i] = Eigen::Matrix<double,1,1>::Zero();
281 zeroSpeedSystem.D[j][i] = Eigen::Matrix<double,1,1>::Zero();
282 }
283 }
284 }
285 RadiationApproximation.insert( std::make_pair(0,zeroSpeedSystem));
286 ZeroSpeedIdentificaitonToCache();
287 }
288
289 void GenerateTimeDomainApproximation()
290 {
291 //std::cout << "Generating time domain state radiation model"<<std::endl;
292 if( IdentificationFromCache() ) return;
293
294 FreqDomainIdent Identification;
295 Identification.SetPotentialData(this);
296 Identification.IdentificationMethod(FreqDomainIdent::VECTOR_FITTING);
297
298 for (unsigned int speed=0; speed < NoSpeed; speed++ )
299 {
300 TimeDomainRadiation speedSystem;
301
302 for( int i=0; i < 6; i++)
303 {
304 for( int j=0; j <= i; j++)
305 {
306 if( (i == 1 || i == 3 || i == 5) && (j==0 || j==2 || j==4) ) continue;
307 if( (j == 1 || j == 3 || j == 5) && (i==0 || i==2 || i==4) ) continue;
308
309 std::cout << "Identifying (i/j) = ("<<speed<<")(" << i << "/" << j << ") ";
310 //CPrecisionTimer T; T.Start();
311 if( is2Ddata() && (i==0||j==0)){
312 FreqDomainIdent::StateSpace SS = Identification.IdentifySS(i,j,speed);
313 speedSystem.Ainf(i,j) = Identification.GetAinf(i,j,speed);
314 speedSystem.A[i][j] = SS.A;
315 speedSystem.B[i][j] = SS.B;
316 speedSystem.C[i][j] = SS.C;
317 speedSystem.D[i][j] = SS.D;
318
319 speedSystem.Ainf(j,i) = speedSystem.Ainf(i,j);
320 speedSystem.A[j][i] = speedSystem.A[i][j];
321 speedSystem.B[j][i] = speedSystem.B[i][j];
322 speedSystem.C[j][i] = speedSystem.C[i][j];
323 speedSystem.D[j][i] = speedSystem.D[i][j];
324 }else{
325 FreqDomainIdent::StateSpace SS = Identification.IdentifySS(i,j);
326 speedSystem.Ainf(i,j) = GetAij(i,j,speed)[0];
327 speedSystem.A[i][j] = Eigen::Matrix<double,1,1>::Zero();
328 speedSystem.B[i][j] = Eigen::Matrix<double,1,1>::Zero();
329 speedSystem.C[i][j] = Eigen::Matrix<double,1,1>::Zero();
330 speedSystem.D[i][j] = Eigen::Matrix<double,1,1>::Zero();
331
332 speedSystem.Ainf(j,i) = GetAij(i,j,speed)[0];
333 speedSystem.A[j][i] = Eigen::Matrix<double,1,1>::Zero();
334 speedSystem.B[j][i] = Eigen::Matrix<double,1,1>::Zero();
335 speedSystem.C[j][i] = Eigen::Matrix<double,1,1>::Zero();
336 speedSystem.D[j][i] = Eigen::Matrix<double,1,1>::Zero();
337 }
338
339 std::cout << "identified new data - " /*<< T.Stop() << " seconds"*/ << std::endl;
340 }
341 }
342 RadiationApproximation.insert( std::make_pair(speed,speedSystem));
343 }
344 IdentificaitonToCache();
345 }
346
347
348 void PrintShip()
349 {
350 Eigen::Matrix<double,6,6> M = GetMassRigidBody();
351
352 std::cout << "%Ship information form file: " << fname << std::endl;
353 std::cout << "%\tLpp:\t" << Lpp << std::endl;
354 std::cout << "%\tB:\t" << Bredth << std::endl;
355 std::cout << "%\tT:\t" << Draught << std::endl;
356 std::cout << "%\tr_cg\t" << LCG<<"\t,\t0\t,\t"<<VCG<< " (from AP and base line)" << std::endl;
357 std::cout << "%\tr_cb\t" << LCB<<"\t,\t0\t,\t"<<VCB<< " (from AP and base line)" << std::endl;
358
359 std::cout << "%\tPotential data set information: " << std::endl;
360 std::cout << "%\t\tNumber of speeds:\t" <<NoSpeed<< std::endl;
361 std::cout << "speed = [\t\t" ;
362 for(unsigned int i=0; i < NoSpeed; i++ ) std::cout << Speed[i] << "\t";
363 std::cout << "];" << std::endl;
364
365 std::cout << "%\t\tNumber of headings:\t" <<NoHeading<< std::endl;
366 std::cout << "head =[\t\t" ;
367 for(unsigned int i=0; i < NoHeading; i++ ) std::cout << Heading[i] << "\t";
368 std::cout << "];" << std::endl;
369
370 std::cout << "%\t\tNumber of frequencies:\t" <<NoFrequency<< std::endl;
371 std::cout << "freq =[\t\t" ;
372 for(unsigned int i=0; i < NoFrequency; i++ ) std::cout << Frequency[i] << "\t";
373 std::cout << "];" << std::endl;
374
375 for(unsigned int i=0; i < NoSpeed; i++ )
376 {
377 std::cout << "%Total mass matrix for Speed " << i << " V="<<Speed[i] << " m/s" << std::endl;
378
379 Eigen::Matrix<double,6,6> Ainf = GetAinf(i);//FIXME: Direction
380
381 Eigen::Matrix<double,6,6> Mtotal = M + Ainf;
382 std::cout << std::scientific;
383 for( int j=0; j < 6; j++)
384 {
385 if( j== 0) std::cout << "M{" <<i+1<<"}\t=[";
386 std::cout <<std::setprecision(4) <<"\t"<< Mtotal(j,0) << "\t," << Mtotal(j,1) << "\t," << Mtotal(j,2) << "\t," << Mtotal(j,3) << "\t," << Mtotal(j,4) << "\t,"<<Mtotal(j,5);
387 if( j== 5) std::cout <<"];";
388 std::cout <<";" <<std::endl;
389 }
390
391 }
392
393
394 for(unsigned int i=0; i < NoSpeed; i++)
395 {
396 std::cout << "%Fluid memory model for Speed " << i << " V="<<Speed[i] << " m/s" << std::endl;
397 TimeDomainRadiation mod = RadiationApproximation[i];
398
399 for( int j=0; j < 6; j++ )
400 {
401 for( int k=0; k < 6; k++ )
402 {
403 size_t num_states = mod.A[k][j].cols();
404
405 if( num_states == 1 )
406 {
407 std::cout << "A{"<<i+1<<"}{"<<j+1<<"}{"<<k+1<<"}" <<"=\t" << "[" << mod.A[j][k](0,0) << "];"<< std::endl;
408 std::cout << std::endl;
409 std::cout << "B{"<<i+1<<"}{"<<j+1<<"}{"<<k+1<<"}" <<" =\t" << "[" << mod.B[j][k](0,0) << "];"<< std::endl;
410 std::cout << std::endl;
411 std::cout << "C{"<<i+1<<"}{"<<j+1<<"}{"<<k+1<<"}" <<" =\t" << "[" << mod.C[j][k](0,0) << "];"<< std::endl;
412 std::cout << std::endl;
413 std::cout << "D{"<<i+1<<"}{"<<j+1<<"}{"<<k+1<<"}" <<" =\t" << "[" << mod.D[j][k](0,0) << "];"<< std::endl;
414 }
415 else if (num_states == 2 )
416 {
417 std::cout << "A{"<<i+1<<"}{"<<j+1<<"}{"<<k+1<<"}"<<" =\t" << "[" << mod.A[j][k](0,0) << "\t,"<<mod.A[j][k](0,1)<<"\t;" << mod.A[j][k](1,0) << "\t,"<<mod.A[j][k](1,1)<<"\t;"<< "];"<< std::endl;
418 std::cout << std::endl;
419 std::cout << "B{"<<i+1<<"}{"<<j+1<<"}{"<<k+1<<"}"<<" =\t" << "[" << mod.B[j][k](0,0) << "\t;" << mod.B[j][k](1,0) <<"];"<< std::endl;
420 std::cout << std::endl;
421 std::cout << "C{"<<i+1<<"}{"<<j+1<<"}{"<<k+1<<"}"<<" =\t" << "[" << mod.C[j][k](0,0) << "\t," << mod.C[j][k](0,1) <<"];"<< std::endl;
422 std::cout << std::endl;
423 std::cout << "D{"<<i+1<<"}{"<<j+1<<"}{"<<k+1<<"}"<<" =\t" << "[" << mod.D[j][k](0,0) <<"];"<< std::endl;
424 }
425 else
426 {
427
428 std::cout << "A{"<<i+1<<"}{"<<j+1<<"}{"<<k+1<<"}"<<" =\t" << "[";
429 for( size_t s=0; s < num_states; s++)
430 {
431 if(s!=0) std::cout <<" \t ";
432 for( size_t c=0; c < num_states; c++)
433 {
434 std::cout << mod.A[j][k](s,c);
435 if( c < (num_states-1) ) std::cout << "\t,";
436 }
437 if( s < (num_states-1) ) std::cout << std::endl;
438 }
439 std::cout << "];"<< std::endl;
440
441 std::cout << std::endl;
442 std::cout << "B{"<<i+1<<"}{"<<j+1<<"}{"<<k+1<<"}"<<" =\t" << "[";
443 for( size_t s=0; s < num_states; s++)
444 {
445 if( s != 0 )std::cout <<" \t ";
446 std::cout << mod.B[j][k](s,0);
447 if( s < (num_states-1) ) std::cout << ";" <<std::endl;
448 }
449 std::cout << "];"<< std::endl;
450 std::cout << std::endl;
451
452 std::cout << "C{"<<i+1<<"}{"<<j+1<<"}{"<<k+1<<"}"<<" =\t" << "[";
453 for( size_t s=0; s < num_states; s++)
454 {
455 std::cout << mod.C[j][k](0,s);
456 if( s < (num_states-1) ) std::cout << ",";
457 }
458 std::cout << "];"<< std::endl;
459
460 std::cout << std::endl;
461 std::cout << "D{"<<i+1<<"}{"<<j+1<<"}{"<<k+1<<"}"<<" =\t" << "[" << mod.D[j][k](0,0) << "];"<< std::endl;
462 }
463 }
464 }
465 }
466 }
467
468 size_t GetTimeDomainStateSize( int DOF1, int DOF2, int SpeedNo=0 )
469 {
470 return RadiationApproximation[SpeedNo].A[DOF1][DOF2].cols();
471 }
472
473 void ZeroSpeedIdentificaitonToCache()
474 {
475 std::string cachebase = fname;
476 std::stringstream ss; ss << cachebase << ".identcache0";
477
478 std::fstream cache(ss.str().c_str(),std::ios::out|std::ios::binary);
479
480 int a_size[2];
481 int b_size[2];
482 int c_size[2];
483 int d_size[2];
484
485 for( int i=0; i < 6; i++)
486 {
487 for( int j=0; j < 6; j++)
488 {
489// a_size[0] = static_cast<int>(RadiationApproximation[0].A[i][j].rows());
490// a_size[1] = static_cast<int>(RadiationApproximation[0].A[i][j].cols());
491
492 a_size[0] = static_cast<int>(RadiationApproximation[0].A[i][j].rows());
493 a_size[1] = static_cast<int>(RadiationApproximation[0].A[i][j].cols());
494
495 cache.write( reinterpret_cast<char*>(a_size), sizeof(a_size) );
496 cache.write( reinterpret_cast<char*>(RadiationApproximation[0].A[i][j].data()),sizeof(double)*a_size[0]*a_size[1] );
497
498// b_size[0] = static_cast<int>(RadiationApproximation[0].B[i][j].rows());
499// b_size[1] = static_cast<int>(RadiationApproximation[0].B[i][j].cols());
500
501 b_size[0] = static_cast<int>(RadiationApproximation[0].B[i][j].rows());
502 b_size[1] = static_cast<int>(RadiationApproximation[0].B[i][j].cols());
503
504 cache.write( reinterpret_cast<char*>(b_size), sizeof(b_size) );
505 cache.write( reinterpret_cast<char*>(RadiationApproximation[0].B[i][j].data()),sizeof(double)*b_size[0]*b_size[1] );
506
507
508 // c_size[0] = static_cast<int>(RadiationApproximation[0].C[i][j].rows());
509 // c_size[1] = static_cast<int>(RadiationApproximation[0].C[i][j].cols());
510
511 c_size[0] = static_cast<int>(RadiationApproximation[0].C[i][j].rows());
512 c_size[1] = static_cast<int>(RadiationApproximation[0].C[i][j].cols());
513
514 cache.write( reinterpret_cast<char*>(c_size), sizeof(c_size) );
515 cache.write( reinterpret_cast<char*>(RadiationApproximation[0].C[i][j].data()),sizeof(double)*c_size[0]*c_size[1] );
516
517// d_size[0] = static_cast<int>(RadiationApproximation[0].D[i][j].rows());
518// d_size[1] = static_cast<int>(RadiationApproximation[0].D[i][j].cols());
519//
520 d_size[0] = static_cast<int>(RadiationApproximation[0].D[i][j].rows());
521 d_size[1] = static_cast<int>(RadiationApproximation[0].D[i][j].cols());
522
523 cache.write( reinterpret_cast<char*>(d_size), sizeof(d_size) );
524 cache.write( reinterpret_cast<char*>(RadiationApproximation[0].D[i][j].data()),sizeof(double)*d_size[0]*d_size[1] );
525
526 //std::cout << "Caching speed: " << MySpeed << " ("<<i<<","<<j<<") - " << " wrote A[ "<<a_size[0]<<" x "<<a_size[1]<<" ]"<<std::endl;
527 }
528 }
529
530 cache.write( reinterpret_cast<char*>(RadiationApproximation[0].Ainf.data()),sizeof(double)*6*6 );
531 cache.close();
532 }
533
534 bool ZeroSpeedIdentificationFromCache()
535 {
536 std::string cachebase = fname;
537 std::stringstream ss; ss << cachebase <<".identcache0";
538 std::fstream cache(ss.str().c_str(),std::ios::in|std::ios::binary);
539 if( !cache.good() ) return false;
540
541 int a_size[2];
542 int b_size[2];
543 int c_size[2];
544 int d_size[2];
545
546 double *A_tmp;
547 double *B_tmp;
548 double *C_tmp;
549 double *D_tmp;
550 double Ainf_tmp[6*6];
551 TimeDomainRadiation radiationModel;
552 for( int i=0; i < 6; i++)
553 {
554 for( int j=0; j < 6; j++)
555 {
556 cache.read( reinterpret_cast<char*>(a_size), sizeof(a_size) );
557 A_tmp = new double[a_size[0]*a_size[1]];
558
559 cache.read( reinterpret_cast<char*>(A_tmp),sizeof(double)*a_size[0]*a_size[1] );
560
561
562 cache.read( reinterpret_cast<char*>(b_size), sizeof(b_size) );
563 B_tmp = new double[b_size[0]*b_size[1]];
564
565 cache.read( reinterpret_cast<char*>(B_tmp),sizeof(double)*b_size[0]*b_size[1] );
566
567 cache.read( reinterpret_cast<char*>(c_size), sizeof(c_size) );
568 C_tmp = new double[c_size[0]*c_size[1]];
569
570 cache.read( reinterpret_cast<char*>(C_tmp),sizeof(double)*c_size[0]*c_size[1] );
571
572
573 cache.read( reinterpret_cast<char*>(d_size), sizeof(d_size) );
574 D_tmp = new double[d_size[0]*d_size[1]];
575
576 cache.read( reinterpret_cast<char*>(D_tmp),sizeof(double)*d_size[0]*d_size[1] );
577
578 //std::cout << "Recovered speed: " << MySpeed << " ("<<i<<","<<j<<") - " << " wrote A["<<a_size[0]<<"x"<<a_size[1]<<"]"<<std::endl;
579
580 radiationModel.A[i][j] = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>::Zero(a_size[0],a_size[1]);
581 radiationModel.B[i][j] = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>::Zero(b_size[0],b_size[1]);
582 radiationModel.C[i][j] = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>::Zero(c_size[0],c_size[1]);
583 radiationModel.D[i][j] = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>::Zero(d_size[0],d_size[1]);
584
585 for( int k=0; k < a_size[0]; k++)
586 for(int l=0; l < a_size[1]; l++)
587 {
588 assert(a_size[0]*a_size[1] > (l+k*a_size[1]) );
589 radiationModel.A[i][j](k,l) = A_tmp[k+l*a_size[0]];
590 }
591 for( int k=0; k < b_size[0]; k++)
592 for(int l=0; l < b_size[1]; l++)
593 {
594 assert(b_size[0]*b_size[1] > (l+k*b_size[1]) );
595 radiationModel.B[i][j](k,l) = B_tmp[k+l*b_size[0]];
596 }
597
598 for( int k=0; k < c_size[0]; k++)
599 for(int l=0; l < c_size[1]; l++)
600 {
601 assert(c_size[0]*c_size[1] > (l+k*c_size[1]) );
602 radiationModel.C[i][j](k,l) = C_tmp[k+l*c_size[0]];
603 }
604 for( int k=0; k < d_size[0]; k++)
605 for(int l=0; l < d_size[1]; l++)
606 {
607 assert(d_size[0]*d_size[1] > (l+k*d_size[1]) );
608 radiationModel.D[i][j](k,l) = D_tmp[k+l*d_size[0]];
609 }
610 delete []A_tmp;
611 delete []B_tmp;
612 delete []C_tmp;
613 delete []D_tmp;
614 }
615 }
616
617 cache.read( reinterpret_cast<char*>(Ainf_tmp),sizeof(double)*6*6);
618
619 for(int i=0; i < 6; i++)
620 for(int j=0; j < 6; j++)
621 radiationModel.Ainf(i,j)=Ainf_tmp[j+i*6];
622
623 radiationModel.Speed = 0;
624 RadiationApproximation.insert( std::make_pair(0,radiationModel) );
625 cache.close();
626 return true;
627 }
628
629 void IdentificaitonToCache()
630 {
631 std::string cachebase = fname;
632 std::stringstream ss; ss << cachebase << ".identcache1";
633
634 std::fstream cache(ss.str().c_str(),std::ios::out|std::ios::binary);
635
636 int a_size[2];
637 int b_size[2];
638 int c_size[2];
639 int d_size[2];
640
641 cache.write( reinterpret_cast<char*>(&NoSpeed),sizeof(unsigned int) );
642
643 for(unsigned int SpeedIdx=0; SpeedIdx < NoSpeed; SpeedIdx++)
644 {
645 double MySpeed = Speed[SpeedIdx];
646
647 cache.write( reinterpret_cast<char*>(&MySpeed),sizeof(double) );
648 for( int i=0; i < 6; i++)
649 {
650 for( int j=0; j < 6; j++)
651 {
652// a_size[0] = static_cast<int>(RadiationApproximation[SpeedIdx].A[i][j].rows());
653// a_size[1] = static_cast<int>(RadiationApproximation[SpeedIdx].A[i][j].cols());
654
655 a_size[0] = static_cast<int>(RadiationApproximation[SpeedIdx].A[i][j].rows());
656 a_size[1] = static_cast<int>(RadiationApproximation[SpeedIdx].A[i][j].cols());
657
658 cache.write( reinterpret_cast<char*>(a_size), sizeof(a_size) );
659 cache.write( reinterpret_cast<char*>(RadiationApproximation[SpeedIdx].A[i][j].data()),sizeof(double)*a_size[0]*a_size[1] );
660
661// b_size[0] = static_cast<int>(RadiationApproximation[SpeedIdx].B[i][j].rows());
662// b_size[1] = static_cast<int>(RadiationApproximation[SpeedIdx].B[i][j].cols());
663
664 b_size[0] = static_cast<int>(RadiationApproximation[SpeedIdx].B[i][j].rows());
665 b_size[1] = static_cast<int>(RadiationApproximation[SpeedIdx].B[i][j].cols());
666
667 cache.write( reinterpret_cast<char*>(b_size), sizeof(b_size) );
668 cache.write( reinterpret_cast<char*>(RadiationApproximation[SpeedIdx].B[i][j].data()),sizeof(double)*b_size[0]*b_size[1] );
669
670// c_size[0] = static_cast<int>(RadiationApproximation[SpeedIdx].C[i][j].rows());
671// c_size[1] = static_cast<int>(RadiationApproximation[SpeedIdx].C[i][j].cols());
672
673 c_size[0] = static_cast<int>(RadiationApproximation[SpeedIdx].C[i][j].rows());
674 c_size[1] = static_cast<int>(RadiationApproximation[SpeedIdx].C[i][j].cols());
675
676 cache.write( reinterpret_cast<char*>(c_size), sizeof(c_size) );
677 cache.write( reinterpret_cast<char*>(RadiationApproximation[SpeedIdx].C[i][j].data()),sizeof(double)*c_size[0]*c_size[1] );
678
679// d_size[0] = static_cast<int>(RadiationApproximation[SpeedIdx].D[i][j].rows());
680// d_size[1] = static_cast<int>(RadiationApproximation[SpeedIdx].D[i][j].cols());
681
682 d_size[0] = static_cast<int>(RadiationApproximation[SpeedIdx].D[i][j].rows());
683 d_size[1] = static_cast<int>(RadiationApproximation[SpeedIdx].D[i][j].cols());
684
685 cache.write( reinterpret_cast<char*>(d_size), sizeof(d_size) );
686 cache.write( reinterpret_cast<char*>(RadiationApproximation[SpeedIdx].D[i][j].data()),sizeof(double)*d_size[0]*d_size[1] );
687
688 //std::cout << "Caching speed: " << MySpeed << " ("<<i<<","<<j<<") - " << " wrote A[ "<<a_size[0]<<" x "<<a_size[1]<<" ]"<<std::endl;
689 }
690 }
691 cache.write( reinterpret_cast<char*>(RadiationApproximation[SpeedIdx].Ainf.data()),sizeof(double)*6*6 );
692 }
693 cache.close();
694 }
695
696 bool IdentificationFromCache()
697 {
698 std::string cachebase = fname;
699 std::stringstream ss; ss << cachebase <<".identcache1";
700 std::fstream cache(ss.str().c_str(),std::ios::in|std::ios::binary);
701 if( !cache.good() ) return false;
702
703 int a_size[2];
704 int b_size[2];
705 int c_size[2];
706 int d_size[2];
707
708 double *A_tmp;
709 double *B_tmp;
710 double *C_tmp;
711 double *D_tmp;
712 double Ainf_tmp[6*6];
713
714 cache.read( reinterpret_cast<char*>(&NoSpeed),sizeof(unsigned int) );
715
716 for(unsigned int SpeedIdx=0; SpeedIdx < NoSpeed; SpeedIdx++)
717 {
718 double MySpeed;
719
720 cache.read( reinterpret_cast<char*>(&MySpeed),sizeof(double) );
721
722 TimeDomainRadiation radiationModel;
723 for( int i=0; i < 6; i++)
724 {
725 for( int j=0; j < 6; j++)
726 {
727 cache.read( reinterpret_cast<char*>(a_size), sizeof(a_size) );
728 A_tmp = new double[a_size[0]*a_size[1]];
729
730 cache.read( reinterpret_cast<char*>(A_tmp),sizeof(double)*a_size[0]*a_size[1] );
731
732 cache.read( reinterpret_cast<char*>(b_size), sizeof(b_size) );
733 B_tmp = new double[b_size[0]*b_size[1]];
734
735 cache.read( reinterpret_cast<char*>(B_tmp),sizeof(double)*b_size[0]*b_size[1] );
736
737 cache.read( reinterpret_cast<char*>(c_size), sizeof(c_size) );
738 C_tmp = new double[c_size[0]*c_size[1]];
739
740 cache.read( reinterpret_cast<char*>(C_tmp),sizeof(double)*c_size[0]*c_size[1] );
741
742 cache.read( reinterpret_cast<char*>(d_size), sizeof(d_size) );
743 D_tmp = new double[d_size[0]*d_size[1]];
744
745 cache.read( reinterpret_cast<char*>(D_tmp),sizeof(double)*d_size[0]*d_size[1] );
746 //std::cout << "Recovered speed: " << MySpeed << " ("<<i<<","<<j<<") - " << " wrote A["<<a_size[0]<<"x"<<a_size[1]<<"]"<<std::endl;
747
748 radiationModel.A[i][j] = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>::Zero(a_size[0],a_size[1]);
749 radiationModel.B[i][j] = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>::Zero(b_size[0],b_size[1]);
750 radiationModel.C[i][j] = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>::Zero(c_size[0],c_size[1]);
751 radiationModel.D[i][j] = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>::Zero(d_size[0],d_size[1]);
752
753 for( int k=0; k < a_size[0]; k++)
754 for(int l=0; l < a_size[1]; l++)
755 {
756 assert(a_size[0]*a_size[1] > (l+k*a_size[1]) );
757 radiationModel.A[i][j](k,l) = A_tmp[k+l*a_size[0]];
758 }
759 for( int k=0; k < b_size[0]; k++)
760 for(int l=0; l < b_size[1]; l++)
761 {
762 assert(b_size[0]*b_size[1] > (l+k*b_size[1]) );
763 radiationModel.B[i][j](k,l) = B_tmp[k+l*b_size[0]];
764 }
765
766 for( int k=0; k < c_size[0]; k++)
767 for(int l=0; l < c_size[1]; l++)
768 {
769 assert(c_size[0]*c_size[1] > (l+k*c_size[1]) );
770 radiationModel.C[i][j](k,l) = C_tmp[k+l*c_size[0]];
771 }
772 for( int k=0; k < d_size[0]; k++)
773 for(int l=0; l < d_size[1]; l++)
774 {
775 assert(d_size[0]*d_size[1] > (l+k*d_size[1]) );
776 radiationModel.D[i][j](k,l) = D_tmp[k+l*d_size[0]];
777 }
778 delete []A_tmp;
779 delete []B_tmp;
780 delete []C_tmp;
781 delete []D_tmp;
782 }
783 }
784
785 cache.read( reinterpret_cast<char*>(Ainf_tmp),sizeof(double)*6*6);
786
787 for(int i=0; i < 6; i++)
788 for(int j=0; j < 6; j++)
789 radiationModel.Ainf(i,j)=Ainf_tmp[j+i*6];
790
791 radiationModel.Speed = MySpeed;
792 RadiationApproximation.insert( std::make_pair(SpeedIdx, radiationModel ));
793 }
794 cache.close();
795 return true;
796 }
797
798
799 void RadiationDerivative(int DOF1, int DOF2, double *derivative, const double *const radiation_state, const double *ship_state, unsigned int SpeedNo=0 )
800 {
801 const TimeDomainRadiation& radiationModel = RadiationApproximation[SpeedNo];
802
803 const size_t state_size = radiationModel.A[DOF1][DOF2].cols();
804 Eigen::Map<const Eigen::MatrixXd,Eigen::Unaligned> X(radiation_state,state_size,1);
805 Eigen::Map<Eigen::MatrixXd,Eigen::Unaligned> dotX(derivative,state_size,1);
806 dotX.noalias() = radiationModel.A[DOF1][DOF2]*X + radiationModel.B[DOF1][DOF2]*ship_state[DOF2];
807 }
808
809 double RadiationOutput(int DOF1, int DOF2, const double *const radiation_state, const double * const ship_state, unsigned int SpeedNo=0)
810 {
811 const TimeDomainRadiation& radiationModel = RadiationApproximation[SpeedNo];
812
813 const size_t state_size = radiationModel.A[DOF1][DOF2].cols();
814 Eigen::Map<const Eigen::MatrixXd,Eigen::Unaligned> X(radiation_state,state_size,1);
815 return (radiationModel.C[DOF1][DOF2]*X + radiationModel.D[DOF1][DOF2]*ship_state[DOF2]).value();
816 }
817
818
819 std::map<int,TimeDomainRadiation,std::less<int>,Eigen::aligned_allocator<std::pair<const int,TimeDomainRadiation> > > RadiationApproximation;
820
821 protected:
822
823 virtual bool is2Ddata(){
824 return false;
825 }
826
827 std::string fname;
828
829 std::vector<std::string> m_warnings;
830
831 // NaN until read from the file
832 double Lpp = std::numeric_limits<double>::quiet_NaN(), Bredth = std::numeric_limits<double>::quiet_NaN(), Draught = std::numeric_limits<double>::quiet_NaN();
833
834 double LCG = std::numeric_limits<double>::quiet_NaN(), VCG = std::numeric_limits<double>::quiet_NaN();
835 double LCB = std::numeric_limits<double>::quiet_NaN(), VCB = std::numeric_limits<double>::quiet_NaN();
836
837
838 unsigned int NoSpeed;
839 unsigned int NoHeading;
840 unsigned int NoDOF;
841 unsigned int NoFrequency;
842
843 double *Speed = nullptr;
844 double *Heading = nullptr;
845 double **Frequency = nullptr;
846
847 FroudeKrylovForce *FroudeFrylov = nullptr;
848 WaveDriftForce *WaveDrift = nullptr;
849
850 Eigen::Matrix<double,6,6> Mrbeig;
851 Eigen::Matrix<double,6,6> **A;
852 Eigen::Matrix<double,6,6> **B;
853 Eigen::Matrix<double,6,6> **C;
854 Eigen::Matrix<double,6,6> **Aadd;
855 Eigen::Matrix<double,6,6> **Badd;
856 Eigen::Matrix<double,6,6> **Cadd;
857
858 Eigen::Matrix<double,3,1> rg;
859 Eigen::Matrix<double,3,1> ro;
860 std::map<std::string,double> param;
861 };
862}
863#endif
Ship::PotentialData
Definition PotentialData.h:29
Ship::PotentialData::Warnings
const std::vector< std::string > & Warnings() const
What the reader noted while it read the data (a missing optional file). The ship logs each once.
Definition PotentialData.h:76
Ship
Simple waypoint object.
Definition CableAttach.h:16
Ship::PotentialData::TimeDomainRadiation
Definition PotentialData.h:34