transformPoints.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 Copyright (C) 2017-2022 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27Application
28 transformPoints
29
30Group
31 grpMeshManipulationUtilities
32
33Description
34 Transforms the mesh points in the polyMesh directory according to the
35 translate, rotate and scale options.
36
37Usage
38 Options are:
39
40 -time value
41 Specify the time to search from and apply the transformation
42 (default is latest)
43
44 -recentre
45 Recentre using the bounding box centre before other operations
46
47 -translate vector
48 Translate the points by the given vector before rotations
49
50 -rotate (vector vector)
51 Rotate the points from the first vector to the second
52
53 -rotate-angle (vector angle)
54 Rotate angle degrees about vector axis.
55
56 -rotate-x angle
57 Rotate (degrees) about x-axis.
58
59 -rotate-y angle
60 Rotate (degrees) about y-axis.
61
62 -rotate-z angle
63 Rotate (degrees) about z-axis.
64
65 or -yawPitchRoll (yawdegrees pitchdegrees rolldegrees)
66 or -rollPitchYaw (rolldegrees pitchdegrees yawdegrees)
67
68 -scale scalar|vector
69 Scale the points by the given scalar or vector on output.
70
71 The any or all of the three options may be specified and are processed
72 in the above order.
73
74 With -rotateFields (in combination with -rotate/yawPitchRoll/rollPitchYaw)
75 it will also read & transform vector & tensor fields.
76
77Note
78 roll (rotation about x)
79 pitch (rotation about y)
80 yaw (rotation about z)
81
82\*---------------------------------------------------------------------------*/
83
84#include "argList.H"
85#include "Time.H"
86#include "fvMesh.H"
87#include "volFields.H"
88#include "surfaceFields.H"
89#include "pointFields.H"
90#include "ReadFields.H"
91#include "regionProperties.H"
92#include "transformField.H"
94#include "axisAngleRotation.H"
96
97using namespace Foam;
98using namespace Foam::coordinateRotations;
99
100// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
101
102template<class GeoField>
103void readAndRotateFields
104(
105 PtrList<GeoField>& flds,
106 const fvMesh& mesh,
107 const dimensionedTensor& rotT,
108 const IOobjectList& objects
109)
110{
111 ReadFields(mesh, objects, flds);
112 for (GeoField& fld : flds)
113 {
114 Info<< "Transforming " << fld.name() << endl;
115 transform(fld, rotT, fld);
116 }
117}
118
119
120void rotateFields
121(
122 const word& regionName,
123 const Time& runTime,
124 const tensor& rotT
125)
126{
127 // Need dimensionedTensor for geometric fields
128 const dimensionedTensor T(rotT);
129
130 #include "createRegionMesh.H"
131
132 // Read objects in time directory
133 IOobjectList objects(mesh, runTime.timeName());
134
135 // Read vol fields.
136
138 readAndRotateFields(vsFlds, mesh, T, objects);
139
141 readAndRotateFields(vvFlds, mesh, T, objects);
142
144 readAndRotateFields(vstFlds, mesh, T, objects);
145
147 readAndRotateFields(vsymtFlds, mesh, T, objects);
148
150 readAndRotateFields(vtFlds, mesh, T, objects);
151
152 // Read surface fields.
153
155 readAndRotateFields(ssFlds, mesh, T, objects);
156
158 readAndRotateFields(svFlds, mesh, T, objects);
159
161 readAndRotateFields(sstFlds, mesh, T, objects);
162
164 readAndRotateFields(ssymtFlds, mesh, T, objects);
165
167 readAndRotateFields(stFlds, mesh, T, objects);
168
169 mesh.write();
170}
171
172
173// Retrieve scaling option
174// - size 0 : no scaling
175// - size 1 : uniform scaling
176// - size 3 : non-uniform scaling
177List<scalar> getScalingOpt(const word& optName, const argList& args)
178{
179 // readListIfPresent handles single or multiple values
180 // - accept 1 or 3 values
181
182 List<scalar> scaling;
183 args.readListIfPresent(optName, scaling);
184
185 if (scaling.size() == 1)
186 {
187 // Uniform scaling
188 }
189 else if (scaling.size() == 3)
190 {
191 // Non-uniform, but may actually be uniform
192 if
193 (
194 equal(scaling[0], scaling[1])
195 && equal(scaling[0], scaling[2])
196 )
197 {
198 scaling.resize(1);
199 }
200 }
201 else if (!scaling.empty())
202 {
204 << "Incorrect number of components, must be 1 or 3." << nl
205 << " -" << optName << ' ' << args[optName].c_str() << endl
206 << exit(FatalError);
207 }
208
209 if (scaling.size() == 1 && equal(scaling[0], 1))
210 {
211 // Scale factor 1 == no scaling
212 scaling.clear();
213 }
214
215 // Zero and negative scaling are permitted
216
217 return scaling;
218}
219
220
221void applyScaling(pointField& points, const List<scalar>& scaling)
222{
223 if (scaling.size() == 1)
224 {
225 Info<< "Scaling points uniformly by " << scaling[0] << nl;
226 points *= scaling[0];
227 }
228 else if (scaling.size() == 3)
229 {
230 Info<< "Scaling points by ("
231 << scaling[0] << ' '
232 << scaling[1] << ' '
233 << scaling[2] << ')' << nl;
234
235 points.replace(vector::X, scaling[0]*points.component(vector::X));
236 points.replace(vector::Y, scaling[1]*points.component(vector::Y));
237 points.replace(vector::Z, scaling[2]*points.component(vector::Z));
238 }
239}
240
241
242// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
243
244int main(int argc, char *argv[])
245{
246 argList::addNote
247 (
248 "Transform (translate / rotate / scale) mesh points.\n"
249 "Note: roll=rotate about x, pitch=rotate about y, yaw=rotate about z"
250 );
251 argList::addOption
252 (
253 "time",
254 "time",
255 "Specify the time to search from and apply the transformation"
256 " (default is latest)"
257 );
258 argList::addBoolOption
259 (
260 "recentre",
261 "Recentre the bounding box before other operations"
262 );
263 argList::addOption
264 (
265 "translate",
266 "vector",
267 "Translate by specified <vector> before rotations"
268 );
269 argList::addBoolOption
270 (
271 "auto-origin",
272 "Use bounding box centre as origin for rotations"
273 );
274 argList::addOption
275 (
276 "origin",
277 "point",
278 "Use specified <point> as origin for rotations"
279 );
280 argList::addOption
281 (
282 "rotate",
283 "(vectorA vectorB)",
284 "Rotate from <vectorA> to <vectorB> - eg, '((1 0 0) (0 0 1))'"
285 );
286 argList::addOption
287 (
288 "rotate-angle",
289 "(vector angle)",
290 "Rotate <angle> degrees about <vector> - eg, '((1 0 0) 45)'"
291 );
292 argList::addOption
293 (
294 "rotate-x", "deg",
295 "Rotate (degrees) about x-axis"
296 );
297 argList::addOption
298 (
299 "rotate-y", "deg",
300 "Rotate (degrees) about y-axis"
301 );
302 argList::addOption
303 (
304 "rotate-z", "deg",
305 "Rotate (degrees) about z-axis"
306 );
307 argList::addOption
308 (
309 "rollPitchYaw",
310 "vector",
311 "Rotate by '(roll pitch yaw)' degrees"
312 );
313 argList::addOption
314 (
315 "yawPitchRoll",
316 "vector",
317 "Rotate by '(yaw pitch roll)' degrees"
318 );
319 argList::addBoolOption
320 (
321 "rotateFields",
322 "Read and transform vector and tensor fields too"
323 );
324 argList::addOption
325 (
326 "scale",
327 "scalar | vector",
328 "Scale by the specified amount - Eg, for uniform [mm] to [m] scaling "
329 "use either '(0.001 0.001 0.001)' or simply '0.001'"
330 );
331
332 // Compatibility with surfaceTransformPoints
333 argList::addOptionCompat("scale", {"write-scale", 0});
334
335 #include "addAllRegionOptions.H"
336 #include "setRootCase.H"
337
338 const bool doRotateFields = args.found("rotateFields");
339
340 // Verify that an operation has been specified
341 {
342 const List<word> operationNames
343 ({
344 "recentre",
345 "translate",
346 "rotate",
347 "rotate-angle",
348 "rotate-x",
349 "rotate-y",
350 "rotate-z",
351 "rollPitchYaw",
352 "yawPitchRoll",
353 "scale"
354 });
355
356 if (!args.count(operationNames))
357 {
359 << "No operation supplied, "
360 << "use at least one of the following:" << nl
361 << " ";
362
363 for (const auto& opName : operationNames)
364 {
366 << " -" << opName;
367 }
368
370 << nl << exit(FatalError);
371 }
372 }
373
374 // ------------------------------------------------------------------------
375
376 #include "createTime.H"
377
378 if (args.found("time"))
379 {
380 if (args["time"] == "constant")
381 {
382 runTime.setTime(instant(0, "constant"), 0);
383 }
384 else
385 {
386 const scalar timeValue = args.get<scalar>("time");
387 runTime.setTime(instant(timeValue), 0);
388 }
389 }
390
391 // Handle -allRegions, -regions, -region
392 #include "getAllRegionOptions.H"
393
394 // ------------------------------------------------------------------------
395
396 forAll(regionNames, regioni)
397 {
398 const word& regionName = regionNames[regioni];
399 const fileName meshDir
400 (
401 polyMesh::regionName(regionName)/polyMesh::meshSubDir
402 );
403
404 if (regionNames.size() > 1)
405 {
406 Info<< "region=" << regionName << nl;
407 }
408
410 (
412 (
413 "points",
414 runTime.findInstance(meshDir, "points"),
415 meshDir,
416 runTime,
417 IOobject::MUST_READ,
418 IOobject::NO_WRITE,
419 false
420 )
421 );
422
423
424 // Begin operations
425
426 vector v;
427 if (args.found("recentre"))
428 {
429 v = boundBox(points).centre();
430 Info<< "Adjust centre " << v << " -> (0 0 0)" << endl;
431 points -= v;
432 }
433
434 if (args.readIfPresent("translate", v))
435 {
436 Info<< "Translating points by " << v << endl;
437 points += v;
438 }
439
440 vector origin;
441 bool useOrigin = args.readIfPresent("origin", origin);
442 if (args.found("auto-origin") && !useOrigin)
443 {
444 useOrigin = true;
445 origin = boundBox(points).centre();
446 }
447
448 if (useOrigin)
449 {
450 Info<< "Set origin for rotations to " << origin << endl;
451 points -= origin;
452 }
453
454
455 // Get a rotation specification
456
457 tensor rot(Zero);
458 bool useRotation(false);
459
460 if (args.found("rotate"))
461 {
462 Pair<vector> n1n2
463 (
464 args.lookup("rotate")()
465 );
466 n1n2[0].normalise();
467 n1n2[1].normalise();
468
469 rot = rotationTensor(n1n2[0], n1n2[1]);
470 useRotation = true;
471 }
472 else if (args.found("rotate-angle"))
473 {
474 const Tuple2<vector, scalar> rotAxisAngle
475 (
476 args.lookup("rotate-angle")()
477 );
478
479 const vector& axis = rotAxisAngle.first();
480 const scalar angle = rotAxisAngle.second();
481
482 Info<< "Rotating points " << nl
483 << " about " << axis << nl
484 << " angle " << angle << nl;
485
486 rot = axisAngle::rotation(axis, angle, true);
487 useRotation = true;
488 }
489 else if (args.found("rotate-x"))
490 {
491 const scalar angle = args.get<scalar>("rotate-x");
492
493 Info<< "Rotating points about x-axis: " << angle << nl;
494
495 rot = axisAngle::rotation(vector::X, angle, true);
496 useRotation = true;
497 }
498 else if (args.found("rotate-y"))
499 {
500 const scalar angle = args.get<scalar>("rotate-y");
501
502 Info<< "Rotating points about y-axis: " << angle << nl;
503
504 rot = axisAngle::rotation(vector::Y, angle, true);
505 useRotation = true;
506 }
507 else if (args.found("rotate-z"))
508 {
509 const scalar angle = args.get<scalar>("rotate-z");
510
511 Info<< "Rotating points about z-axis: " << angle << nl;
512
513 rot = axisAngle::rotation(vector::Z, angle, true);
514 useRotation = true;
515 }
516 else if (args.readIfPresent("rollPitchYaw", v))
517 {
518 Info<< "Rotating points by" << nl
519 << " roll " << v.x() << nl
520 << " pitch " << v.y() << nl
521 << " yaw " << v.z() << nl;
522
523 rot = euler::rotation(euler::eulerOrder::ROLL_PITCH_YAW, v, true);
524 useRotation = true;
525 }
526 else if (args.readIfPresent("yawPitchRoll", v))
527 {
528 Info<< "Rotating points by" << nl
529 << " yaw " << v.x() << nl
530 << " pitch " << v.y() << nl
531 << " roll " << v.z() << nl;
532
533 rot = euler::rotation(euler::eulerOrder::YAW_PITCH_ROLL, v, true);
534 useRotation = true;
535 }
536
537 if (useRotation)
538 {
539 Info<< "Rotating points by " << rot << endl;
540 transform(points, rot, points);
541
542 if (doRotateFields)
543 {
544 rotateFields(regionName, runTime, rot);
545 }
546 }
547
548 // Output scaling
549 applyScaling(points, getScalingOpt("scale", args));
550
551 if (useOrigin)
552 {
553 Info<< "Unset origin for rotations from " << origin << endl;
554 points += origin;
555 }
556
557
558 // Set the precision of the points data to 10
559 IOstream::defaultPrecision(max(10u, IOstream::defaultPrecision()));
560
561 Info<< "Writing points into directory "
562 << runTime.relativePath(points.path()) << nl
563 << endl;
564 points.write();
565 }
566
567 Info<< nl << "End" << nl << endl;
568
569 return 0;
570}
571
572
573// ************************************************************************* //
Field reading functions for post-processing utilities.
Y[inertIndex] max(0.0)
Info<< nl<< "Wrote faMesh in vtk format: "<< writer.output().name()<< nl;}{ vtk::lineWriter writer(aMesh.points(), aMesh.edges(), fileName(aMesh.mesh().time().globalPath()/"finiteArea-edges"));writer.writeGeometry();writer.beginCellData(4);writer.writeProcIDs();{ Field< scalar > fld(faMeshTools::flattenEdgeField(aMesh.magLe(), true))
A primitive field of type <T> with automated input and output.
Definition: IOField.H:58
List of IOobjects with searching and retrieving facilities.
Definition: IOobjectList.H:59
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:170
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: List.H:77
void resize(const label len)
Adjust allocated size of list.
Definition: ListI.H:139
void clear()
Clear the list, i.e. set size to zero.
Definition: ListI.H:116
An ordered pair of two objects of type <T> with first() and second() elements.
Definition: Pair.H:69
A list of pointers to objects of type <T>, with allocation/deallocation management of the pointers....
Definition: PtrList.H:73
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
Definition: Time.H:80
A 2-tuple for storing two objects of dissimilar types. The container is similar in purpose to std::pa...
Definition: Tuple2.H:58
bool empty() const noexcept
True if the UList is empty (ie, size() is zero)
Definition: UListI.H:427
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
Extract command arguments and options from the supplied argc and argv parameters.
Definition: argList.H:124
T get(const label index) const
Get a value from the argument at index.
Definition: argListI.H:278
bool readListIfPresent(const word &optName, List< T > &list) const
Definition: argListI.H:394
label count(const UList< word > &optionNames) const
Return how many of the specified options were used.
Definition: argList.C:1760
bool found(const word &optName) const
Return true if the named option is found.
Definition: argListI.H:178
bool readIfPresent(const word &optName, T &val) const
Read a value from the named option if present.
Definition: argListI.H:323
ITstream lookup(const word &optName) const
Return an input stream from the named option.
Definition: argListI.H:184
A bounding box defined in terms of min/max extrema points.
Definition: boundBox.H:64
point centre() const
The centre (midpoint) of the bounding box.
Definition: boundBoxI.H:115
Generic dimensioned Type class.
A class for handling file names.
Definition: fileName.H:76
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:91
An instant of time. Contains the time value and name. Uses Foam::Time when formatting the name.
Definition: instant.H:56
Tensor of scalars, i.e. Tensor<scalar>.
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
A class for handling words, derived from Foam::string.
Definition: word.H:68
const volScalarField & T
dynamicFvMesh & mesh
engineTime & runTime
Foam::word regionName(Foam::polyMesh::defaultRegion)
Required Variables.
wordList regionNames
const pointField & points
bool equal(const Matrix< Form1, Type > &A, const Matrix< Form2, Type > &B, const bool verbose=false, const label maxDiffs=10, const scalar relTol=1e-5, const scalar absTol=1e-8)
Compare matrix elements for absolute or relative equality.
Definition: MatrixTools.C:34
Namespace for coordinate system rotations.
Definition: axesRotation.C:38
Namespace for OpenFOAM.
wordList ReadFields(const typename GeoMesh::Mesh &mesh, const IOobjectList &objects, PtrList< GeometricField< Type, PatchField, GeoMesh > > &fields, const bool syncPar=true, const bool readOldTime=false)
Read Geometric fields of templated type.
dimensionSet transform(const dimensionSet &ds)
Return the argument; transformations do not change the dimensions.
Definition: dimensionSet.C:536
tensor rotationTensor(const vector &n1, const vector &n2)
Rotational transformation tensor from vector n1 to n2.
Definition: transform.H:51
messageStream Info
Information stream (stdout output on master, null elsewhere)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
error FatalError
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
Foam::argList args(argc, argv)
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
Foam::surfaceFields.
Spatial transformation functions for primitive fields.
Spatial transformation functions for GeometricField.