CppNoddy  0.92
Loading...
Searching...
No Matches
Utility.h
Go to the documentation of this file.
1/// \file Utility.h
2/// A spec for a collection of utility functions.
3
4#ifndef UTILITY_H
5#define UTILITY_H
6
7#include <OneD_Node_Mesh.h>
8#include <string>
9#include <ctime>
10#include <numeric>
11#include <algorithm>
12#include <numeric>
13#include <sys/types.h>
14#include <sys/stat.h>
15
16#include <Types.h>
17#include <Exceptions.h>
18#include <FortranBLAS.h>
19#include <TwoD_Node_Mesh.h>
20#include <TwoD_Mapped_Node_Mesh.h>
21#include <Sequential_Matrix_base.h>
22
23namespace CppNoddy {
24 /// Some utility methods associated with CppNoddy containers.
25 namespace Utility {
26
27 /// Return a DENSE vector with the nodal points of a uniform
28 /// mesh distributed between the upper/lower bounds as specified
29 /// \param lower The lower bound of the uniform nodal distribution
30 /// \param upper The upper bound of the uniform nodal distribution
31 /// \param N The number of nodal points
32 DenseVector<double> uniform_node_vector(const double& lower, const double& upper, const std::size_t& N);
33
34 /// Return a DENSE vector with the nodal points of a non-uniform
35 /// mesh distributed between the upper/lower bounds as specified
36 /// with more nodes clustered near lower or upper depending upon
37 /// the differencee of the power from unity. When power=1 this should
38 /// provide a uniform mesh.
39 /// \param lower The lower bound of the uniform nodal distribution
40 /// \param upper The upper bound of the uniform nodal distribution
41 /// \param N The number of nodal points
42 /// \param power A measure of the non-uniformity
43 /// \return A vector of nodal positions with a power law distribution
44 DenseVector<double> power_node_vector(const double& lower, const double& upper, const std::size_t& N, const double& power);
45
46 /// Return a dense vector with two uniform distributions in two separate
47 /// regions.
48 /// \param lower The first node
49 /// \param mid The node that defines the boundary between the uniform meshes
50 /// \param upper The final node
51 /// \param N1 The number of nodes in the first region
52 /// \param N2 The number of nodes in the second region
53 /// \return A combined vector of nodes of length N1+N2
54 DenseVector<double> two_uniform_node_vector(const double& lower, const double& mid, const double& upper, const std::size_t& N1, const std::size_t& N2);
55
56 /// Return a dense vector with two uniform distributions in two separate
57 /// regions.
58 /// \param lower The first node
59 /// \param mid1 The node that defines the first interior boundary
60 /// \param mid2 The node that defines the second interior boundary
61 /// \param upper The final node
62 /// \param N1 The number of nodes in the first region
63 /// \param N2 The number of nodes in the second region
64 /// \param N3 The number of nodes in the third region
65 /// \return A combined vector of nodes of length N1+N2+N3
66 DenseVector<double> three_uniform_node_vector(const double& lower, const double& mid1, const double& mid2, const double& upper, const std::size_t& N1, const std::size_t& N2, const std::size_t& N3);
67
68 /// Return a dense vector of nodal positions with more nodes concentrated
69 /// at the mid point of the range.
70 /// \param lower The first nodal position.
71 /// \param upper The final nodal position.
72 /// \param N The number of nodes required.
73 /// \param power A measure of the non-uniformity, power = 1 => uniform distribution
74 DenseVector<double> mid_weighted_node_vector(const double& lower, const double& upper, const std::size_t& N, const double& power);
75
76 //
77 //
78 // SOME typical ops
79 //
80 //
81
82 template <typename _Type>
83 void vels_from_streamfn_Cartesian(const TwoD_Node_Mesh<_Type>& source, TwoD_Node_Mesh<_Type>& uv) {
84 std::cout << "PLAIN MESH version\n";
85 std::size_t Nx(source.get_nnodes().first);
86 std::size_t Ny(source.get_nnodes().second);
87 double dx(source.coord(1,1).first - source.coord(0,0).first);
88 double dy(source.coord(1,1).second - source.coord(0,0).second);
89 // differentiate the streamfunction to get the velocity field
90 {
91 // west internal nodes
92 std::size_t i(0);
93 for(std::size_t j = 1; j < Ny - 1; ++j) {
94 uv(i, j, 0) = (source(i, j + 1, 0) - source(i, j - 1, 0)) / (2 * dy);
95 uv(i, j, 1) = -(-source(i + 2, j, 0) + 4. * source(i + 1, j, 0) - 3. * source(i, j, 0)) / (2 * dx);
96 }
97 }
98 {
99 // east internal nodes
100 std::size_t i(Nx - 1);
101 for(std::size_t j = 1; j < Ny - 1; ++j) {
102 uv(i, j, 0) = (source(i, j + 1, 0) - source(i, j - 1, 0)) / (2 * dy);
103 uv(i, j, 1) = -(source(i - 2, j, 0) - 4. * source(i - 1, j, 0) + 3. * source(i, j, 0)) / (2 * dx);
104 }
105 }
106 {
107 // south internal nodes
108 std::size_t j(0);
109 for(std::size_t i = 1; i < Nx - 1; ++i) {
110 uv(i, j, 0) = (-source(i, j + 2, 0) + 4. * source(i, j + 1, 0) - 3. * source(i, j, 0)) / (2 * dy);
111 uv(i, j, 1) = -(source(i + 1, j, 0) - source(i - 1, j, 0)) / (2 * dx);
112 }
113 }
114 {
115 // north internal nodes
116 std::size_t j(Ny - 1);
117 for(std::size_t i = 1; i < Nx - 1; ++i) {
118 uv(i, j, 0) = (source(i, j - 2, 0) - 4. * source(i, j - 1, 0) + 3. * source(i, j, 0)) / (2 * dy);
119 uv(i, j, 1) = -(source(i + 1, j, 0) - source(i - 1, j, 0)) / (2 * dx);
120 }
121 }
122 {
123 // corner nodes
124 {
125 // sw
126 std::size_t i(0);
127 std::size_t j(0);
128 uv(i, j, 0) = (-source(i, j + 2, 0) + 4. * source(i, j + 1, 0) - 3. * source(i, j, 0)) / (2 * dy);
129 uv(i, j, 1) = -(-source(i + 2, j, 0) + 4. * source(i + 1, j, 0) - 3. * source(i, j, 0)) / (2 * dx);
130 }
131 {
132 // nw
133 std::size_t i(0);
134 std::size_t j(Ny - 1);
135 uv(i, j, 0) = (source(i, j - 2, 0) - 4. * source(i, j - 1, 0) + 3. * source(i, j, 0)) / (2 * dy);
136 uv(i, j, 1) = -(-source(i + 2, j, 0) + 4. * source(i + 1, j, 0) - 3. * source(i, j, 0)) / (2 * dx);
137 }
138 {
139 // ne
140 std::size_t i(Nx - 1);
141 std::size_t j(Ny - 1);
142 uv(i, j, 0) = (source(i, j - 2, 0) - 4. * source(i, j - 1, 0) + 3. * source(i, j, 0)) / (2 * dy);
143 uv(i, j, 1) = -(source(i - 2, j, 0) - 4. * source(i - 1, j, 0) + 3. * source(i, j, 0)) / (2 * dx);
144 }
145 {
146 // se
147 std::size_t i(Nx - 1);
148 std::size_t j(0);
149 uv(i, j, 0) = (-source(i, j + 2, 0) + 4. * source(i, j + 1, 0) - 3. * source(i, j, 0)) / (2 * dy);
150 uv(i, j, 1) = -(source(i - 2, j, 0) - 4. * source(i - 1, j, 0) + 3. * source(i, j, 0)) / (2 * dx);
151 }
152 }
153 {
154 // interior nodes
155 for(std::size_t i = 1; i < Nx - 1; ++i) {
156 for(std::size_t j = 1; j < Ny - 1; ++j) {
157 uv(i, j, 0) = (source(i, j + 1, 0) - source(i, j - 1, 0)) / (2 * dy);
158 uv(i, j, 1) = -(source(i + 1, j, 0) - source(i - 1, j, 0)) / (2 * dx);
159 }
160 }
161 }
162 }
163
164
165 template <typename _Type>
166 void vels_from_streamfn_Cartesian(const TwoD_Mapped_Node_Mesh<_Type>& source, TwoD_Mapped_Node_Mesh<_Type>& uv) {
167 std::cout << "MAPPED MESH version\n";
168 std::size_t Nx(source.get_nnodes().first);
169 std::size_t Ny(source.get_nnodes().second);
170 double dX( source.get_comp_step_sizes().first );
171 double dY( source.get_comp_step_sizes().second );
172 // differentiate the streamfunction to get the velocity field
173 {
174 // west internal nodes
175 std::size_t i(0);
176 for(std::size_t j = 1; j < Ny - 1; ++j) {
177 double x = source.coord(i,j).first;
178 double Xd = source.FnComp_Xd(x);
179 double y = source.coord(i,j).second;
180 double Yd = source.FnComp_Yd(y);
181 uv(i, j, 0) = (source(i, j + 1, 0) - source(i, j - 1, 0)) * Yd / (2 * dY);
182 uv(i, j, 1) = -(-source(i + 2, j, 0) + 4. * source(i + 1, j, 0) - 3. * source(i, j, 0)) * Xd / (2 * dX);
183 }
184 }
185 {
186 // east internal nodes
187 std::size_t i(Nx - 1);
188 for(std::size_t j = 1; j < Ny - 1; ++j) {
189 double x = source.coord(i,j).first;
190 double Xd = source.FnComp_Xd(x);
191 double y = source.coord(i,j).second;
192 double Yd = source.FnComp_Yd(y);
193 uv(i, j, 0) = (source(i, j + 1, 0) - source(i, j - 1, 0)) * Yd / (2 * dY);
194 uv(i, j, 1) = -(source(i - 2, j, 0) - 4. * source(i - 1, j, 0) + 3. * source(i, j, 0)) * Xd / (2 * dX);
195 }
196 }
197 {
198 // south internal nodes
199 std::size_t j(0);
200 for(std::size_t i = 1; i < Nx - 1; ++i) {
201 double x = source.coord(i,j).first;
202 double Xd = source.FnComp_Xd(x);
203 double y = source.coord(i,j).second;
204 double Yd = source.FnComp_Yd(y);
205 uv(i, j, 0) = (-source(i, j + 2, 0) + 4. * source(i, j + 1, 0) - 3. * source(i, j, 0))* Yd / (2 * dY);
206 uv(i, j, 1) = -(source(i + 1, j, 0) - source(i - 1, j, 0)) * Xd / (2 * dX);
207 }
208 }
209 {
210 // north internal nodes
211 std::size_t j(Ny - 1);
212 for(std::size_t i = 1; i < Nx - 1; ++i) {
213 double x = source.coord(i,j).first;
214 double Xd = source.FnComp_Xd(x);
215 double y = source.coord(i,j).second;
216 double Yd = source.FnComp_Yd(y);
217 uv(i, j, 0) = (source(i, j - 2, 0) - 4. * source(i, j - 1, 0) + 3. * source(i, j, 0)) * Yd / (2 * dY);
218 uv(i, j, 1) = -(source(i + 1, j, 0) - source(i - 1, j, 0)) * Xd / (2 * dX);
219 }
220 }
221 {
222 // corner nodes
223 {
224 // sw
225 std::size_t i(0);
226 std::size_t j(0);
227 double x = source.coord(i,j).first;
228 double Xd = source.FnComp_Xd(x);
229 double y = source.coord(i,j).second;
230 double Yd = source.FnComp_Yd(y);
231 uv(i, j, 0) = (-source(i, j + 2, 0) + 4. * source(i, j + 1, 0) - 3. * source(i, j, 0))* Xd / (2 * dY);
232 uv(i, j, 1) = -(-source(i + 2, j, 0) + 4. * source(i + 1, j, 0) - 3. * source(i, j, 0))* Yd / (2 * dY);
233 }
234 {
235 // nw
236 std::size_t i(0);
237 std::size_t j(Ny - 1);
238 double x = source.coord(i,j).first;
239 double Xd = source.FnComp_Xd(x);
240 double y = source.coord(i,j).second;
241 double Yd = source.FnComp_Yd(y);
242 uv(i, j, 0) = (source(i, j - 2, 0) - 4. * source(i, j - 1, 0) + 3. * source(i, j, 0)) * Yd / (2 * dY);
243 uv(i, j, 1) = -(-source(i + 2, j, 0) + 4. * source(i + 1, j, 0) - 3. * source(i, j, 0))* Xd / (2 * dX);
244 }
245 {
246 // ne
247 std::size_t i(Nx - 1);
248 std::size_t j(Ny - 1);
249 double x = source.coord(i,j).first;
250 double Xd = source.FnComp_Xd(x);
251 double y = source.coord(i,j).second;
252 double Yd = source.FnComp_Yd(y);
253 uv(i, j, 0) = (source(i, j - 2, 0) - 4. * source(i, j - 1, 0) + 3. * source(i, j, 0)) * Yd / (2 * dY);
254 uv(i, j, 1) = -(source(i - 2, j, 0) - 4. * source(i - 1, j, 0) + 3. * source(i, j, 0)) * Xd/ (2 * dX);
255 }
256 {
257 // se
258 std::size_t i(Nx - 1);
259 std::size_t j(0);
260 double x = source.coord(i,j).first;
261 double Xd = source.FnComp_Xd(x);
262 double y = source.coord(i,j).second;
263 double Yd = source.FnComp_Yd(y);
264 uv(i, j, 0) = (-source(i, j + 2, 0) + 4. * source(i, j + 1, 0) - 3. * source(i, j, 0))* Yd / (2 * dY);
265 uv(i, j, 1) = -(source(i - 2, j, 0) - 4. * source(i - 1, j, 0) + 3. * source(i, j, 0)) * Xd / (2 * dX);
266 }
267 }
268 {
269 // interior nodes
270 for(std::size_t i = 1; i < Nx - 1; ++i) {
271 for(std::size_t j = 1; j < Ny - 1; ++j) {
272 double x = source.coord(i,j).first;
273 double Xd = source.FnComp_Xd(x);
274 double y = source.coord(i,j).second;
275 double Yd = source.FnComp_Yd(y);
276 uv(i, j, 0) = (source(i, j + 1, 0) - source(i, j - 1, 0)) * Yd / (2 * dY);
277 uv(i, j, 1) = -(source(i + 1, j, 0) - source(i - 1, j, 0)) * Xd / (2 * dX);
278 }
279 }
280 }
281 }
282
283
284
285
286 //
287 //
288 // MATRIX OPERATIONS
289 //
290 //
291
292 /// BLAS wrapper to do DOUBLE DENSE A_{MxK} * B_{KxN} = C_{MxN}
293 /// Since this is a Fortran library, it assumes a column_major
294 /// format, but CppNoddy uses row_major. To be consistent we'll
295 /// simply do (B^T)_{NxK} * (A^T)_{KxM} = (C^T)_{NxM} instead.
296 /// Note that inversion of the transpose of the result C^T
297 /// is handled implicitly via the construction of C.
298 /// \param A First dense double matrix to be multiplied
299 /// \param B Second dense double matrix to be multiplied
300 /// \return The result of the multiplication C=A*B
301 DenseMatrix<double> multiply(DenseMatrix<double>& A, DenseMatrix<double>& B);
302
303 //
304 //
305 // VECTOR OPERATIONS
306 //
307 //
308
309 /// Templated dot product.
310 /// \param X First dense vector
311 /// \param Y Second dense vector
312 /// \return The dot product
313 template <typename _Type>
314 _Type dot(const DenseVector<_Type>& X, const DenseVector<_Type>& Y) {
315 if(X.size() != Y.size()) {
316 std::string problem;
317 problem = "The Utilities::dot method has been called \n";
318 problem += "with two unequal length vectors.";
319 throw ExceptionGeom(problem, X.size(), Y.size());
320 }
321 return inner_product(X.begin(), X.end(), Y.begin(), _Type(0.0));
322 }
323
324 template <typename _Type>
325 int sgn(const _Type& a) {
326 if(a > (_Type)0) {
327 return 1;
328 } else if(a < (_Type)0) {
329 return -1;
330 } else {
331 return 0;
332 }
333 }
334
335 //
336 //
337 // COMPLEX UTILS
338 //
339 //
340
341 /// Return a double DENSE vector containing the real part
342 /// of a complex DENSE vector
343 /// \param X The complex vector to take the real part of
344 DenseVector<double> real(const DenseVector<D_complex>& X);
345
346 /// Return a double DENSE vector containing the imaginary part
347 /// of a complex DENSE vector
348 /// \param X The complex vector to take the imaginary part of
349 DenseVector<double> imag(const DenseVector<D_complex>& X);
350
351 //
352 //
353 // STRING TWEAKERY
354 //
355 //
356
357 /// Return an integer value as a string - useful for file naming
358 /// \param val The integer value to be stringified.
359 std::string stringify(const int &val);
360
361 /// Return a double value as a string - useful for file naming.
362 /// \param val The double value to be stringified
363 /// \param p Precision to be used in the output
364 std::string stringify(const double &val, int p);
365
366
367
368 double max_abs_location( OneD_Node_Mesh<double> &mesh , unsigned var);
369
370 double max_abs_location_range( OneD_Node_Mesh<double> &mesh , unsigned var, double left, double right);
371
372 }
373
374} // end namespace
375
376#endif // UTILITY_H
Exceptions.h
The collection of CppNoddy exceptions.
FortranBLAS.h
Some interface definitions for calling external fortran routines from BLAS.
OneD_Node_Mesh.h
A specification for a one dimensional mesh object.
Sequential_Matrix_base.h
A base matrix class to ensure a consistent interface between the inheriting dense/banded matrix class...
TwoD_Mapped_Node_Mesh.h
A specification for a two dimensional mapped mesh object.
TwoD_Node_Mesh.h
A specification for a two dimensional mesh object.
Types.h
Some standard typedefs.
CppNoddy::DenseMatrix
A matrix class that constructs a DENSE matrix as a row major std::vector of DenseVectors.
Definition: DenseMatrix.h:25
CppNoddy::DenseVector
An DenseVector class – a dense vector object.
Definition: DenseVector.h:34
CppNoddy::DenseVector::size
std::size_t size() const
A pass-thru definition to get the size of the vector.
Definition: DenseVector.h:330
CppNoddy::DenseVector::begin
elt_iter begin()
Pass through to the storage container.
Definition: DenseVector.h:82
CppNoddy::DenseVector::end
elt_iter end()
Pass through to the storage container.
Definition: DenseVector.h:102
CppNoddy::ExceptionGeom
An exception class to be thrown when a container of incorrect geometry used in any class/method.
Definition: Exceptions.h:47
CppNoddy::OneD_Node_Mesh
A one dimensional mesh utility object.
Definition: OneD_Node_Mesh.h:24
CppNoddy::TwoD_Mapped_Node_Mesh
A two dimensional (mapped) mesh utility object.
Definition: TwoD_Mapped_Node_Mesh.h:21
CppNoddy::TwoD_Node_Mesh
A two dimensional mesh utility object.
Definition: TwoD_Node_Mesh.h:19
CppNoddy::Utility::real
DenseVector< double > real(const DenseVector< D_complex > &X)
Return a double DENSE vector containing the real part of a complex DENSE vector.
Definition: Utility.cpp:177
CppNoddy::Utility::multiply
DenseMatrix< double > multiply(DenseMatrix< double > &A, DenseMatrix< double > &B)
BLAS wrapper to do DOUBLE DENSE A_{MxK} * B_{KxN} = C_{MxN} Since this is a Fortran library,...
Definition: Utility.cpp:225
CppNoddy::Utility::dot
_Type dot(const DenseVector< _Type > &X, const DenseVector< _Type > &Y)
Templated dot product.
Definition: Utility.h:314
CppNoddy::Utility::sgn
int sgn(const _Type &a)
Definition: Utility.h:325
CppNoddy::Utility::three_uniform_node_vector
DenseVector< double > three_uniform_node_vector(const double &lower, const double &mid1, const double &mid2, const double &upper, const std::size_t &N1, const std::size_t &N2, const std::size_t &N3)
Return a dense vector with two uniform distributions in two separate regions.
Definition: Utility.cpp:141
CppNoddy::Utility::power_node_vector
DenseVector< double > power_node_vector(const double &lower, const double &upper, const std::size_t &N, const double &power)
Return a DENSE vector with the nodal points of a non-uniform mesh distributed between the upper/lower...
Definition: Utility.cpp:123
CppNoddy::Utility::mid_weighted_node_vector
DenseVector< double > mid_weighted_node_vector(const double &lower, const double &upper, const std::size_t &N, const double &power)
Return a dense vector of nodal positions with more nodes concentrated at the mid point of the range.
Definition: Utility.cpp:155
CppNoddy::Utility::vels_from_streamfn_Cartesian
void vels_from_streamfn_Cartesian(const TwoD_Node_Mesh< _Type > &source, TwoD_Node_Mesh< _Type > &uv)
Definition: Utility.h:83
CppNoddy::Utility::two_uniform_node_vector
DenseVector< double > two_uniform_node_vector(const double &lower, const double &mid, const double &upper, const std::size_t &N1, const std::size_t &N2)
Return a dense vector with two uniform distributions in two separate regions.
Definition: Utility.cpp:131
CppNoddy::Utility::max_abs_location
double max_abs_location(OneD_Node_Mesh< double > &mesh, unsigned var)
Definition: Utility.cpp:54
CppNoddy::Utility::imag
DenseVector< double > imag(const DenseVector< D_complex > &X)
Return a double DENSE vector containing the imaginary part of a complex DENSE vector.
Definition: Utility.cpp:185
CppNoddy::Utility::max_abs_location_range
double max_abs_location_range(OneD_Node_Mesh< double > &mesh, unsigned var, double left, double right)
Definition: Utility.cpp:82
CppNoddy::Utility::stringify
std::string stringify(const int &val)
Return an integer value as a string - useful for file naming.
Definition: Utility.cpp:193
CppNoddy::Utility::uniform_node_vector
DenseVector< double > uniform_node_vector(const double &lower, const double &upper, const std::size_t &N)
Return a DENSE vector with the nodal points of a uniform mesh distributed between the upper/lower bou...
Definition: Utility.cpp:113
CppNoddy
A collection of OO numerical routines aimed at simple (typical) applied problems in continuum mechani...
Definition: ArcLength_base.h:13
Generated on Tue Jan 1 1980 00:00:00 for CppNoddy by   doxygen 1.9.5

© 2012

R.E. Hewitt