Skip to content

Commit bafc14f

Browse files
ebertolazziebertolazzi
committed
corrected reading of SplineSurf and build of BiCubic and BiQuintic
1 parent 2fa22ac commit bafc14f

19 files changed

+314
-358
lines changed

‎.travis.yml‎

Lines changed: 0 additions & 62 deletions
This file was deleted.

‎src/SplineAkima.cc‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -135,7 +135,7 @@ namespace Splines {
135135

136136
do {
137137
// cerca intervallo monotono strettamente crescente
138-
while ( ++iend < m_npts && m_X[iend-1] < m_X[iend] ) {}
138+
for ( ++iend; iend < m_npts && m_X[iend-1] < m_X[iend]; ++iend ) {}
139139
Akima_build( m_X+ibegin, m_Y+ibegin, m_Yp+ibegin, m_work, iend-ibegin );
140140
ibegin = iend;
141141
} while ( iend < m_npts );

‎src/SplineAkima2D.cc‎

Lines changed: 21 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -156,23 +156,29 @@ namespace Splines {
156156
Akima2Dspline::write_to_stream( ostream_type & s ) const {
157157
fmt::print( s, "Nx={} Ny={}\n", m_nx, m_ny );
158158
for ( integer i{1}; i < m_nx; ++i ) {
159+
real_type const dx{ m_X[i] - m_X[i-1] };
159160
for ( integer j{1}; j < m_ny; ++j ) {
160-
integer const i00 { ipos_C(i-1,j-1) };
161-
integer const i10 { ipos_C(i,j-1) };
162-
integer const i01 { ipos_C(i-1,j) };
163-
integer const i11 { ipos_C(i,j) };
161+
integer const i00 { ipos_C(i-1,j-1) };
162+
integer const i10 { ipos_C(i,j-1) };
163+
integer const i01 { ipos_C(i-1,j) };
164+
integer const i11 { ipos_C(i,j) };
165+
real_type const dy { m_Y[j] - m_Y[j-1] };
164166
fmt::print( s,
165-
"patch ({},{)\n"
166-
" DX = {:<12.4} DY = {:<12.4}\n"
167-
" Z00 = {:<12.4} Z01 = {:<12.4} Z10 = {:<12.4} Z11 = {:<12.4}\n"
168-
" Dx00 = {:<12.4} Dx01 = {:<12.4} Dx10 = {:<12.4} Dx11 = {:<12.4}\n"
169-
" Dy00 = {:<12.4} Dy01 = {:<12.4} Dy10 = {:<12.4} Dy11 = {:<12.4}\n",
170-
i, j,
171-
m_X[i] - m_X[i-1],
172-
m_Y[j] - m_Y[j-1],
173-
m_Z[i00], m_Z[i01], m_Z[i10], m_Z[i11],
174-
m_DX[i00], m_DX[i01], m_DX[i10], m_DX[i11],
175-
m_DY[i00], m_DY[i01], m_DY[i10], m_DY[i11]
167+
"patch ({},{})\n"
168+
" DX = {:<12.4} DY = {:<12.4}\n"
169+
" Z00 = {:<12.4} Z10 = {:<12.4}\n"
170+
" Z01 = {:<12.4} Z11 = {:<12.4}\n"
171+
" Dx00 = {:<12.4} Dx10 = {:<12.4}\n"
172+
" Dx01 = {:<12.4} Dx11 = {:<12.4}\n"
173+
" Dy00 = {:<12.4} Dy10 = {:<12.4}\n"
174+
" Dy01 = {:<12.4} Dy11 = {:<12.4}\n"
175+
" Dxy00 = {:<12.4} Dxy10 = {:<12.4}\n"
176+
" Dxy01 = {:<12.4} Dxy11 = {:<12.4}\n",
177+
i, j, dx, dy,
178+
m_Z[i00], m_Z[i10], m_Z[i01], m_Z[i11],
179+
m_DX[i00], m_DX[i10], m_DX[i01], m_DX[i11],
180+
m_DY[i00], m_DY[i10], m_DY[i01], m_DY[i11],
181+
m_DXY[i00], m_DXY[i10], m_DXY[i01], m_DXY[i11]
176182
);
177183
}
178184
}

‎src/SplineBessel.cc‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -109,7 +109,7 @@ namespace Splines {
109109
integer iend{0};
110110
do {
111111
// cerca intervallo monotono strettamente crescente
112-
while ( ++iend < m_npts && m_X[iend-1] < m_X[iend] ) {}
112+
for ( ++iend; iend < m_npts && m_X[iend-1] < m_X[iend]; ++iend ) {}
113113
Bessel_build( m_X+ibegin, m_Y+ibegin, m_Yp+ibegin, iend-ibegin );
114114
ibegin = iend;
115115
} while ( iend < m_npts );

‎src/SplineBiCubic.cc‎

Lines changed: 20 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -85,23 +85,29 @@ namespace Splines {
8585
BiCubicSpline::write_to_stream( ostream_type & s ) const {
8686
fmt::print( "Nx = {} Ny = {}\n", m_nx, m_ny );
8787
for ( integer i{1}; i < m_nx; ++i ) {
88+
real_type dx{ m_X[i]-m_X[i-1] };
8889
for ( integer j{1}; j < m_ny; ++j ) {
89-
integer const i00 { this->ipos_C(i-1,j-1) };
90-
integer const i10 { this->ipos_C(i,j-1) };
91-
integer const i01 { this->ipos_C(i-1,j) };
92-
integer const i11 { this->ipos_C(i,j) };
90+
integer const i00 { ipos_C(i-1,j-1) };
91+
integer const i10 { ipos_C(i,j-1) };
92+
integer const i01 { ipos_C(i-1,j) };
93+
integer const i11 { ipos_C(i,j) };
94+
real_type const dy { m_Y[j]-m_Y[j-1] };
9395
fmt::print( s,
9496
"patch ({},{})\n"
95-
"DX = {:<12.4} DY = {:<12.4}\n"
96-
"Z00 = {:<12.4} Z01 = {:<12.4} Z10 = {:<12.4} Z11 = {:<12.4}\n"
97-
"Dx00 = {:<12.4} Dx01 = {:<12.4} Dx10 = {:<12.4} Dx11 = {:<12.4}\n"
98-
"Dy00 = {:<12.4} Dy01 = {:<12.4} Dy10 = {:<12.4} Dy11 = {:<12.4}\n",
99-
i, j,
100-
m_X[i]-m_X[i-1],
101-
m_Y[j]-m_Y[j-1],
102-
m_Z[i00], m_Z[i01], m_Z[i10], m_Z[i11],
103-
m_DX[i00], m_DX[i01], m_DX[i10], m_DX[i11],
104-
m_DY[i00], m_DY[i01], m_DY[i10], m_DY[i11]
97+
" DX = {:<12.4} DY = {:<12.4}\n"
98+
" Z00 = {:<12.4} Z10 = {:<12.4}\n"
99+
" Z01 = {:<12.4} Z11 = {:<12.4}\n"
100+
" Dx00 = {:<12.4} Dx10 = {:<12.4}\n"
101+
" Dx01 = {:<12.4} Dx11 = {:<12.4}\n"
102+
" Dy00 = {:<12.4} Dy10 = {:<12.4}\n"
103+
" Dy01 = {:<12.4} Dy11 = {:<12.4}\n"
104+
" Dxy00 = {:<12.4} Dxy10 = {:<12.4}\n"
105+
" Dxy01 = {:<12.4} Dxy11 = {:<12.4}\n",
106+
i, j, dx, dy,
107+
m_Z[i00], m_Z[i10], m_Z[i01], m_Z[i11],
108+
m_DX[i00], m_DX[i10], m_DX[i01], m_DX[i11],
109+
m_DY[i00], m_DY[i10], m_DY[i01], m_DY[i11],
110+
m_DXY[i00], m_DXY[i10], m_DXY[i01], m_DXY[i11]
105111
);
106112
}
107113
}

‎src/SplineBiQuintic.cc‎

Lines changed: 77 additions & 63 deletions
Original file line numberDiff line numberDiff line change
@@ -40,63 +40,71 @@ namespace Splines {
4040
BiQuinticSpline::make_spline() {
4141

4242
integer const dim{ m_nx*m_ny };
43-
mem.reallocate( 8*dim );
44-
m_DX = mem( dim );
45-
m_DY = mem( dim );
46-
m_DXY = mem( dim );
47-
m_DXX = mem( dim );
48-
m_DYY = mem( dim );
49-
m_DXYY = mem( dim );
50-
m_DXXY = mem( dim );
51-
m_DXXYY = mem( dim );
52-
53-
// calcolo derivate
54-
QuinticSpline sp, sp1;
55-
//
43+
m_mem_biquintic.reallocate( 8*dim );
44+
m_DX = m_mem_biquintic( dim );
45+
m_DY = m_mem_biquintic( dim );
46+
m_DXY = m_mem_biquintic( dim );
47+
m_DXX = m_mem_biquintic( dim );
48+
m_DYY = m_mem_biquintic( dim );
49+
m_DXYY = m_mem_biquintic( dim );
50+
m_DXXY = m_mem_biquintic( dim );
51+
m_DXXYY = m_mem_biquintic( dim );
52+
53+
auto minmod = [] ( real_type a, real_type b ) -> real_type {
54+
if ( a*b <= 0 ) return 0;
55+
if ( a > 0 ) return std::min(a,b);
56+
return std::max(a,b);
57+
};
58+
59+
PchipSpline sp;
60+
61+
// calcolo derivate DX, DY, DXY
5662
for ( integer j{0}; j < m_ny; ++j ) {
57-
sp.build( m_X, 1, &m_Z[ this->ipos_C(0,j) ], m_ny, m_nx );
58-
for ( integer i{0}; i < m_nx; ++i ) {
59-
integer const ij{ this->ipos_C(i,j) };
60-
m_DX[ij] = sp.yp_node(i);
61-
m_DXX[ij] = sp.ypp_node(i);
62-
}
63+
sp.build( m_X, 1, &z_node_ref(0,j), m_ny, m_nx );
64+
for ( integer i{0}; i < m_nx; ++i ) Dx_node_ref(i,j) = sp.yp_node(i);
6365
}
6466
for ( integer i{0}; i < m_nx; ++i ) {
65-
sp.build( m_Y, 1, &m_Z[ this->ipos_C(i,0) ], 1, m_ny );
66-
for ( integer j{0}; j < m_ny; ++j ) {
67-
integer const ij{ this->ipos_C(i,j) };
68-
m_DY[ij] = sp.yp_node(j);
69-
m_DYY[ij] = sp.ypp_node(j);
70-
}
67+
sp.build( m_Y, 1, &z_node_ref(i,0), 1, m_ny );
68+
for ( integer j{0}; j < m_ny; ++j ) Dy_node_ref(i,j) = sp.yp_node(j);
69+
}
70+
71+
// mixed
72+
for ( integer j{0}; j < m_ny; ++j ) {
73+
sp.build( m_X, 1, &Dx_node_ref(0,j), m_ny, m_nx );
74+
for ( integer i{0}; i < m_nx; ++i ) Dxy_node_ref(i,j) = sp.yp_node(i);
7175
}
72-
// interpolate derivative
7376
for ( integer i{0}; i < m_nx; ++i ) {
74-
sp.build( m_Y, 1, &m_DX[ this->ipos_C(i,0) ], 1, m_ny );
75-
sp1.build( m_Y, 1, &m_DXX[ this->ipos_C(i,0) ], 1, m_ny );
76-
for ( integer j{0}; j < m_ny; ++j ) {
77-
integer const ij{ this->ipos_C(i,j) };
78-
m_DXY[ij] = sp.yp_node(j);
79-
m_DXYY[ij] = sp.ypp_node(j);
80-
m_DXXY[ij] = sp1.yp_node(j);
81-
m_DXXYY[ij] = sp1.ypp_node(j);
82-
}
77+
sp.build( m_Y, 1, &Dy_node_ref(i,0), 1, m_ny );
78+
for ( integer j{0}; j < m_ny; ++j ) Dxy_node_ref(i,j) = minmod( Dxy_node_ref(i,j), sp.yp_node(j) );
8379
}
84-
// interpolate derivative again
80+
81+
// calcolo derivate DXX, DYY, DXXY, DXYY, DXXYY
8582
for ( integer j{0}; j < m_ny; ++j ) {
86-
sp.build( m_X, 1, &m_DY[ this->ipos_C(0,j) ], m_ny, m_nx );
87-
sp1.build( m_X, 1, &m_DYY[ this->ipos_C(0,j) ], m_ny, m_nx );
88-
for ( integer i{0}; i < m_nx; ++i ) {
89-
integer const ij{ this->ipos_C(i,j) };
90-
m_DXY[ij] += sp.yp_node(i); m_DXY[ij] /= 2;
91-
m_DXXY[ij] += sp.ypp_node(i); m_DXXY[ij] /= 2;
92-
m_DXYY[ij] += sp1.yp_node(i); m_DXYY[ij] /= 2;
93-
m_DXXYY[ij] += sp1.ypp_node(i); m_DXXYY[ij] /= 2;
94-
}
83+
sp.build( m_X, 1, &Dx_node_ref(0,j), m_ny, m_nx );
84+
for ( integer i{0}; i < m_nx; ++i ) Dxx_node_ref(i,j) = sp.yp_node(i);
85+
}
86+
for ( integer i{0}; i < m_nx; ++i ) {
87+
sp.build( m_Y, 1, &Dy_node_ref(i,0), 1, m_ny );
88+
for ( integer j{0}; j < m_ny; ++j ) Dyy_node_ref(i,j) = sp.yp_node(j);
89+
}
90+
for ( integer j{0}; j < m_ny; ++j ) {
91+
sp.build( m_X, 1, &Dxy_node_ref(0,j), m_ny, m_nx );
92+
for ( integer i{0}; i < m_nx; ++i ) Dxxy_node_ref(i,j) = sp.yp_node(i);
93+
}
94+
for ( integer i{0}; i < m_nx; ++i ) {
95+
sp.build( m_Y, 1, &Dxy_node_ref(i,0), 1, m_ny );
96+
for ( integer j{0}; j < m_ny; ++j ) Dxyy_node_ref(i,j) = sp.yp_node(j);
97+
}
98+
// mixed
99+
for ( integer j{0}; j < m_ny; ++j ) {
100+
sp.build( m_X, 1, &Dxyy_node_ref(0,j), m_ny, m_nx );
101+
for ( integer i{0}; i < m_nx; ++i ) Dxxyy_node_ref(i,j) = sp.yp_node(i);
102+
}
103+
for ( integer i{0}; i < m_nx; ++i ) {
104+
sp.build( m_Y, 1, &Dxxy_node_ref(i,0), 1, m_ny );
105+
for ( integer j{0}; j < m_ny; ++j ) Dxxy_node_ref(i,j) = minmod( Dxxy_node_ref(i,j), sp.yp_node(j) );
95106
}
96107

97-
//std::fill( DXY.begin(), DXY.end(), 0 );
98-
//std::fill( DXX.begin(), DXX.end(), 0 );
99-
//std::fill( DYY.begin(), DYY.end(), 0 );
100108
}
101109

102110
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
@@ -105,23 +113,29 @@ namespace Splines {
105113
BiQuinticSpline::write_to_stream( ostream_type & s ) const {
106114
fmt::print( s, "Nx = {} Ny = {}\n", m_nx, m_ny );
107115
for ( integer i{1}; i < m_nx; ++i ) {
116+
real_type dx{ m_X[i]-m_X[i-1] };
108117
for ( integer j{1}; j < m_ny; ++j ) {
109-
integer const i00 { this->ipos_C(i-1,j-1) };
110-
integer const i10 { this->ipos_C(i,j-1) };
111-
integer const i01 { this->ipos_C(i-1,j) };
112-
integer const i11 { this->ipos_C(i,j) };
118+
integer const i00 { ipos_C(i-1,j-1) };
119+
integer const i10 { ipos_C(i,j-1) };
120+
integer const i01 { ipos_C(i-1,j) };
121+
integer const i11 { ipos_C(i,j) };
122+
real_type const dy { m_Y[j]-m_Y[j-1] };
113123
fmt::print( s,
114-
"patch({},{})\n"
115-
" DX = {:<12.4} DY = {:<12.4}\n"
116-
" Z00 = {:<12.4} Z01 = {:<12.4} Z10 = {:<12.4} Z11 = {:<12.4}\n"
117-
" Dx00 = {:<12.4} Dx01 = {:<12.4} Dx10 = {:<12.4} Dx10 = {:<12.4}\n"
118-
" Dy00 = {:<12.4} Dy01 = {:<12.4} Dy10 = {:<12.4} Dy11 = {:<12.4}\n",
119-
i, j,
120-
m_X[i]-m_X[i-1],
121-
m_Y[j]-m_Y[j-1],
122-
m_Z[i00], m_Z[i01], m_Z[i10], m_Z[i11],
123-
m_DX[i00], m_DX[i01], m_DX[i10], m_DX[i11],
124-
m_DY[i00], m_DY[i01], m_DY[i10], m_DY[i11]
124+
"patch ({},{})\n"
125+
" DX = {:<12.4} DY = {:<12.4}\n"
126+
" Z00 = {:<12.4} Z10 = {:<12.4}\n"
127+
" Z01 = {:<12.4} Z11 = {:<12.4}\n"
128+
" Dx00 = {:<12.4} Dx10 = {:<12.4}\n"
129+
" Dx01 = {:<12.4} Dx11 = {:<12.4}\n"
130+
" Dy00 = {:<12.4} Dy10 = {:<12.4}\n"
131+
" Dy01 = {:<12.4} Dy11 = {:<12.4}\n"
132+
" Dxy00 = {:<12.4} Dxy10 = {:<12.4}\n"
133+
" Dxy01 = {:<12.4} Dxy11 = {:<12.4}\n",
134+
i, j, dx, dy,
135+
m_Z[i00], m_Z[i10], m_Z[i01], m_Z[i11],
136+
m_DX[i00], m_DX[i10], m_DX[i01], m_DX[i11],
137+
m_DY[i00], m_DY[i10], m_DY[i01], m_DY[i11],
138+
m_DXY[i00], m_DXY[i10], m_DXY[i01], m_DXY[i11]
125139
);
126140
}
127141
}

‎src/SplineBilinear.cc‎

Lines changed: 8 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -123,18 +123,19 @@ namespace Splines {
123123
fmt::print( s, "Nx={} Ny={}\n", m_nx, m_ny );
124124
for ( integer i{1}; i < m_nx; ++i ) {
125125
for ( integer j{1}; j < m_ny; ++j ) {
126-
integer const i00 { this->ipos_C(i-1,j-1) };
127-
integer const i10 { this->ipos_C(i,j-1) };
128-
integer const i01 { this->ipos_C(i-1,j) };
129-
integer const i11 { this->ipos_C(i,j) };
126+
integer const i00 { ipos_C(i-1,j-1) };
127+
integer const i10 { ipos_C(i,j-1) };
128+
integer const i01 { ipos_C(i-1,j) };
129+
integer const i11 { ipos_C(i,j) };
130130
fmt::print( s,
131131
"patch ({},{})\n"
132-
" DX = {:<12.4} DY = {:<12.4}\n"
133-
" Z00 = {:<12.4} Z01 = {:<12.4} Z10 = {:<12.4} Z11 = {:<12.4}\n",
132+
" DX = {:<12.4} DY = {:<12.4}\n"
133+
" Z00 = {:<12.4} Z10 = {:<12.4}\n"
134+
" Z01 = {:<12.4} Z11 = {:<12.4}\n",
134135
i, j,
135136
m_X[i]-m_X[i-1],
136137
m_Y[j]-m_Y[j-1],
137-
m_Z[i00], m_Z[i01], m_Z[i10], m_Z[i11]
138+
m_Z[i00], m_Z[i10], m_Z[i01], m_Z[i11]
138139
);
139140
}
140141
}

0 commit comments

Comments
 (0)