CatmullRomSpline.C
Go to the documentation of this file.
1/*---------------------------------------------------------------------------*\
2 ========= |
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4 \\ / O peration |
5 \\ / A nd | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2011-2016 OpenFOAM Foundation
9-------------------------------------------------------------------------------
10License
11 This file is part of OpenFOAM.
12
13 OpenFOAM is free software: you can redistribute it and/or modify it
14 under the terms of the GNU General Public License as published by
15 the Free Software Foundation, either version 3 of the License, or
16 (at your option) any later version.
17
18 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
19 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21 for more details.
22
23 You should have received a copy of the GNU General Public License
24 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
25
26\*---------------------------------------------------------------------------*/
27
28#include "CatmullRomSpline.H"
29
30// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
31
32Foam::scalar Foam::CatmullRomSpline::derivative
33(
34 const label segment,
35 const scalar mu
36) const
37{
38 const point& p0 = points()[segment];
39 const point& p1 = points()[segment+1];
40
41 // determine the end points
42 point e0;
43 point e1;
44
45 if (segment == 0)
46 {
47 // end: simple reflection
48 e0 = 2*p0 - p1;
49 }
50 else
51 {
52 e0 = points()[segment-1];
53 }
54
55 if (segment+1 == nSegments())
56 {
57 // end: simple reflection
58 e1 = 2*p1 - p0;
59 }
60 else
61 {
62 e1 = points()[segment+2];
63 }
64 const point derivativePoint
65 (
66 0.5 *
67 (
68 (-e0 + p1)
69 + mu *
70 (
71 2 * (2*e0 - 5*p0 + 4*p1 - e1)
72 + mu * 3 * (-e0 + 3*p0 - 3*p1 + e1)
73 )
74 )
75 );
76 return mag(derivativePoint);
77}
78
79
80// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
81
83(
84 const pointField& knots,
85 const bool closed
86)
87:
88 polyLine(knots, closed)
89{}
90
91
92// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
93
95{
96 // endpoints
97 if (mu < SMALL)
98 {
99 return points().first();
100 }
101 else if (mu > 1 - SMALL)
102 {
103 return points().last();
104 }
105
106 scalar lambda = mu;
107 label segment = localParameter(lambda);
108 return position(segment, lambda);
109}
110
111
113(
114 const label segment,
115 const scalar mu
116) const
117{
118 // out-of-bounds
119 if (segment < 0)
120 {
121 return points().first();
122 }
123 else if (segment > nSegments())
124 {
125 return points().last();
126 }
127
128 const point& p0 = points()[segment];
129 const point& p1 = points()[segment+1];
130
131 // special cases - no calculation needed
132 if (mu <= 0.0)
133 {
134 return p0;
135 }
136 else if (mu >= 1.0)
137 {
138 return p1;
139 }
140
141
142 // determine the end points
143 point e0;
144 point e1;
145
146 if (segment == 0)
147 {
148 // end: simple reflection
149 e0 = 2*p0 - p1;
150 }
151 else
152 {
153 e0 = points()[segment-1];
154 }
155
156 if (segment+1 == nSegments())
157 {
158 // end: simple reflection
159 e1 = 2*p1 - p0;
160 }
161 else
162 {
163 e1 = points()[segment+2];
164 }
165
166
167 return
168 0.5 *
169 (
170 (2*p0)
171 + mu *
172 (
173 (-e0 + p1)
174 + mu *
175 (
176 (2*e0 - 5*p0 + 4*p1 - e1)
177 + mu*(-e0 + 3*p0 - 3*p1 + e1)
178 )
179 )
180 );
181}
182
183
185{
186 const solveScalar xi[5]=
187 {
188 -0.9061798459386639927976,
189 -0.5384693101056830910363,
190 0,
191 0.5384693101056830910363,
192 0.9061798459386639927976
193 };
194 const solveScalar wi[5]=
195 {
196 0.2369268850561890875143,
197 0.4786286704993664680413,
198 0.5688888888888888888889,
199 0.4786286704993664680413,
200 0.2369268850561890875143
201 };
202 scalar sum=0;
203 for (label segment=0;segment<nSegments();segment++)
204 {
205 for (int i=0;i<5;i++)
206 {
207 sum+=wi[i]*derivative(segment,(xi[i]+1.0)/2.0)/2.0;
208 }
209 }
210
211 return sum;
212}
213
214
215// ************************************************************************* //
An implementation of Catmull-Rom splines (sometimes known as Overhauser splines).
scalar length() const
The length of the curve.
An ordered pair of two objects of type <T> with first() and second() elements.
Definition: Pair.H:69
T & first()
Return the first element of the list.
Definition: UListI.H:202
T & last()
Return the last element of the list.
Definition: UListI.H:216
A series of straight line segments, which can also be interpreted as a series of control points for s...
Definition: polyLine.H:56
const pointField & points() const noexcept
Return const-access to the control-points.
Definition: polyLine.C:112
label nSegments() const noexcept
The number of line segments.
Definition: polyLine.C:118
const volScalarField & mu
const volScalarField & p0
Definition: EEqn.H:36
const pointField & points
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
vector point
Point is a vector.
Definition: point.H:43
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
Pair< point > segment
dimensionedScalar lambda("lambda", dimTime/sqr(dimLength), laminarTransport)