Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F88407181
cubic_spline.hh
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Fri, Oct 18, 15:50
Size
7 KB
Mime Type
text/x-c++
Expires
Sun, Oct 20, 15:50 (2 d)
Engine
blob
Format
Raw Data
Handle
21766744
Attached To
rLIBMULTISCALE LibMultiScale
cubic_spline.hh
View Options
/**
* @file cubic_spline.hh
*
* @author Srinivasa Babu Ramisetti <srinivasa.ramisetti@epfl.ch>
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @date Wed Dec 11 12:36:48 2013
*
* @brief This is a cubic spline interpolator
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* LibMultiScale is free software: you can redistribute it and/or modify it
* under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* LibMultiScale is distributed in the hope that it will be useful, but
* WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with LibMultiScale. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef __LIBMULTISCALE_CUBIC_SPLINE_HH__
#define __LIBMULTISCALE_CUBIC_SPLINE_HH__
/* -------------------------------------------------------------------------- */
#include <cmath>
#include <vector>
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
class
CubicSpline
{
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public
:
CubicSpline
()
{
bc_first_point
=
0.0
;
bc_last_point
=
0.0
;
spacing
=
1.0
;
spacing_inverse
=
1.0
;
sixinv
=
1.0
/
6.0
;
};
~
CubicSpline
(){};
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public
:
//! build the cubic spline coefficients
void
buildSplineCoefficients
(
std
::
vector
<
Real
>
&
x
,
std
::
vector
<
Real
>
&
f
,
Real
spacing
);
//! evaluate the function at x using cubic spline interoplation
void
evaluateSpline
(
Real
x
,
Real
&
y
);
//! fast evaluate the function at x using cubic spline interoplation works
//! only for equally spaced data points
void
fastEvaluateSpline
(
Real
x
,
Real
&
y
);
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public
:
//! set boundary conditions at first and last points
void
setBCs
(
Real
x0
,
Real
xn
);
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
// private:
protected
:
//! vector to store the x values
std
::
vector
<
Real
>
x_f
;
//! vector to store the f(x) values
std
::
vector
<
Real
>
poly_f
;
//! vector to store the coefficient b
std
::
vector
<
Real
>
b
;
//! vector to store the coefficient c
std
::
vector
<
Real
>
c
;
//! vector to store the coefficient d
std
::
vector
<
Real
>
d
;
//! vector to store the second derivative of f(x) i.e. df_dyy
std
::
vector
<
Real
>
df_dyy
;
//! boundar condition at first data point
Real
bc_first_point
;
//! boundar condition at last data point
Real
bc_last_point
;
//! equal spacing between points
Real
spacing
;
//! equal inverse spacing between points
Real
spacing_inverse
;
//! variable to store constant factor
Real
sixinv
;
//! number of data points
UInt
n
;
};
/* -------------------------------------------------------------------------- */
/* inline functions */
/* -------------------------------------------------------------------------- */
//! build the cubic spline coefficients for not a knot boundary conditions
inline
void
CubicSpline
::
buildSplineCoefficients
(
std
::
vector
<
Real
>
&
x
,
std
::
vector
<
Real
>
&
f
,
Real
spacing
)
{
Real
p
,
qn
,
sig
,
un
;
std
::
vector
<
Real
>
u
;
n
=
x
.
size
();
// internal variables for later use
this
->
spacing
=
spacing
;
spacing_inverse
=
1.
/
spacing
;
x_f
=
x
;
poly_f
=
f
;
u
.
assign
(
n
-
1
,
0.0
);
df_dyy
.
assign
(
n
,
0.0
);
c
.
assign
(
n
,
0.0
);
b
.
assign
(
n
,
0.0
);
d
.
assign
(
n
,
0.0
);
// set the boundary condition at first point
df_dyy
[
0
]
=
-
0.5
;
u
[
0
]
=
(
3.0
/
(
x
[
1
]
-
x
[
0
]))
*
((
f
[
1
]
-
f
[
0
])
/
(
x
[
1
]
-
x
[
0
])
-
bc_first_point
);
for
(
UInt
i
=
1
;
i
<
n
-
1
;
i
++
)
{
sig
=
(
x
[
i
]
-
x
[
i
-
1
])
/
(
x
[
i
+
1
]
-
x
[
i
-
1
]);
p
=
sig
*
df_dyy
[
i
-
1
]
+
2.0
;
df_dyy
[
i
]
=
(
sig
-
1.0
)
/
p
;
u
[
i
]
=
(
f
[
i
+
1
]
-
f
[
i
])
/
(
x
[
i
+
1
]
-
x
[
i
])
-
(
f
[
i
]
-
f
[
i
-
1
])
/
(
x
[
i
]
-
x
[
i
-
1
]);
u
[
i
]
=
(
6.0
*
u
[
i
]
/
(
x
[
i
+
1
]
-
x
[
i
-
1
])
-
sig
*
u
[
i
-
1
])
/
p
;
}
// set the boundary condition at last point
qn
=
0.5
;
un
=
(
3.0
/
(
x
[
n
-
1
]
-
x
[
n
-
2
]))
*
(
bc_last_point
-
(
f
[
n
-
1
]
-
f
[
n
-
2
])
/
(
x
[
n
-
1
]
-
x
[
n
-
2
]));
// store second derivatives of f in df_dyy
df_dyy
[
n
-
1
]
=
(
un
-
qn
*
u
[
n
-
2
])
/
(
qn
*
df_dyy
[
n
-
2
]
+
1.0
);
for
(
int
k
=
n
-
2
;
k
>=
0
;
k
--
)
{
df_dyy
[
k
]
=
df_dyy
[
k
]
*
df_dyy
[
k
+
1
]
+
u
[
k
];
}
// compute cubic spline polynomial coeffs b,c,d
for
(
UInt
i
=
0
;
i
<
n
;
++
i
)
c
[
i
]
=
0.5
*
df_dyy
[
i
];
for
(
UInt
i
=
0
;
i
<
n
-
1
;
++
i
)
{
Real
h
=
(
x
[
i
+
1
]
-
x
[
i
]);
b
[
i
]
=
(
f
[
i
+
1
]
-
f
[
i
])
/
h
-
h
*
(
c
[
i
+
1
]
+
2
*
c
[
i
])
/
3.0
;
d
[
i
]
=
(
c
[
i
+
1
]
-
c
[
i
])
/
(
h
*
3.0
);
}
}
//! evaluate the function at x using cubic spline interoplation (useful for not
//! equally spaced points)
inline
void
CubicSpline
::
evaluateSpline
(
Real
x
,
Real
&
y
)
{
static
UInt
pklo
=
0
,
pkhi
=
1
;
Real
h
,
b
,
a
;
UInt
klo
,
khi
,
k
;
if
(
x_f
[
pklo
]
<=
x
&&
x_f
[
pkhi
]
>
x
)
{
klo
=
pklo
;
khi
=
pkhi
;
}
else
{
klo
=
0
;
khi
=
n
-
1
;
while
(
khi
-
klo
>
1
)
{
k
=
(
khi
+
klo
)
>>
1
;
if
(
x_f
[
k
]
>
x
)
khi
=
k
;
else
klo
=
k
;
}
}
h
=
x_f
[
khi
]
-
x_f
[
klo
];
if
(
h
==
0
)
{
LM_FATAL
(
"This should not happen. Cubic spline cannot be evaluated at this "
"point "
<<
x
);
}
// compute the interpolated value
a
=
(
x_f
[
khi
]
-
x
)
/
h
;
b
=
(
x
-
x_f
[
klo
])
/
h
;
y
=
a
*
poly_f
[
klo
]
+
b
*
poly_f
[
khi
]
+
((
a
*
a
*
a
-
a
)
*
df_dyy
[
klo
]
+
(
b
*
b
*
b
-
b
)
*
df_dyy
[
khi
])
*
(
h
*
h
)
*
sixinv
;
}
inline
void
CubicSpline
::
fastEvaluateSpline
(
Real
x
,
Real
&
y
)
{
UInt
i
=
0
;
// LM_ASSERT(x < x_f[0],"Cubic spline interpolation cannot be evaluated for
// an external point\n");
i
=
(
x
-
x_f
[
0
])
*
spacing_inverse
;
if
(
i
>=
x_f
.
size
()
-
1
)
{
y
=
0.0
;
return
;
}
Real
h
=
x
-
x_f
[
i
];
// compute the interpolated value;
y
=
poly_f
[
i
]
+
h
*
(
b
[
i
]
+
h
*
(
c
[
i
]
+
h
*
d
[
i
]));
LM_ASSERT
(
!
std
::
isnan
(
y
),
"problem"
);
LM_ASSERT
(
!
std
::
isinf
(
y
),
"problem"
);
}
inline
void
CubicSpline
::
setBCs
(
Real
x0
,
Real
xn
)
{
bc_first_point
=
x0
;
bc_last_point
=
xn
;
}
__END_LIBMULTISCALE__
#endif
/* __LIBMULTISCALE_CUBIC_SPLINE_HH__ */
Event Timeline
Log In to Comment