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