surfaceLambdaMuSmooth.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) 2015-2021 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 surfaceLambdaMuSmooth
29
30Group
31 grpSurfaceUtilities
32
33Description
34 Smooth a surface using lambda/mu smoothing.
35
36 To get laplacian smoothing, set lambda to the relaxation factor and mu to
37 zero.
38
39 Provide an edgeMesh file containing points that are not to be moved during
40 smoothing in order to preserve features.
41
42 lambda/mu smoothing: G. Taubin, IBM Research report Rc-19923 (02/01/95)
43 "A signal processing approach to fair surface design"
44
45\*---------------------------------------------------------------------------*/
46
47#include "argList.H"
48#include "boundBox.H"
49#include "edgeMesh.H"
50#include "matchPoints.H"
51#include "MeshedSurfaces.H"
52
53using namespace Foam;
54
55// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
56
58(
59 const meshedSurface& s,
60 const bitSet& fixedPoints
61)
62{
63 const labelListList& pointEdges = s.pointEdges();
64
65 tmp<pointField> tavg(new pointField(s.nPoints(), Zero));
66 pointField& avg = tavg.ref();
67
68 forAll(pointEdges, vertI)
69 {
70 vector& avgPos = avg[vertI];
71
72 if (fixedPoints.test(vertI))
73 {
74 avgPos = s.localPoints()[vertI];
75 }
76 else
77 {
78 const labelList& pEdges = pointEdges[vertI];
79
80 forAll(pEdges, myEdgeI)
81 {
82 const edge& e = s.edges()[pEdges[myEdgeI]];
83
84 label otherVertI = e.otherVertex(vertI);
85
86 avgPos += s.localPoints()[otherVertI];
87 }
88
89 avgPos /= pEdges.size();
90 }
91 }
92
93 return tavg;
94}
95
96
97void getFixedPoints
98(
99 const edgeMesh& feMesh,
100 const pointField& points,
101 bitSet& fixedPoints
102)
103{
104 scalarList matchDistance(feMesh.points().size(), 1e-1);
105 labelList from0To1;
106
107 bool matchedAll = matchPoints
108 (
109 feMesh.points(),
110 points,
111 matchDistance,
112 false,
113 from0To1
114 );
115
116 if (!matchedAll)
117 {
119 << "Did not match all feature points to points on the surface"
120 << endl;
121 }
122
123 forAll(from0To1, fpI)
124 {
125 if (from0To1[fpI] != -1)
126 {
127 fixedPoints.set(from0To1[fpI]);
128 }
129 }
130}
131
132
133// Main program:
134
135int main(int argc, char *argv[])
136{
137 argList::addNote
138 (
139 "Smooth a surface using lambda/mu smoothing.\n"
140 "For laplacian smoothing, set lambda to the relaxation factor"
141 " and mu to zero."
142 );
143
144 argList::noParallel();
145 argList::validOptions.clear();
146 argList::addArgument("input", "The input surface file");
147 argList::addArgument("lambda", "On the interval [0,1]");
148 argList::addArgument("mu", "On the interval [0,1]");
149 argList::addArgument("iterations", "The number of iterations to perform");
150 argList::addArgument("output", "The output surface file");
151
152 argList::addOption
153 (
154 "featureFile",
155 "Fix points from a file containing feature points and edges"
156 );
157 argList args(argc, argv);
158
159 const auto surfFileName = args.get<fileName>(1);
160 const auto lambda = args.get<scalar>(2);
161 const auto mu = args.get<scalar>(3);
162 const auto iters = args.get<label>(4);
163 const auto outFileName = args.get<fileName>(5);
164
165 if (lambda < 0 || lambda > 1)
166 {
168 << lambda << endl
169 << "0: no change 1: move vertices to average of neighbours"
170 << exit(FatalError);
171 }
172 if (mu < 0 || mu > 1)
173 {
175 << mu << endl
176 << "0: no change 1: move vertices to average of neighbours"
177 << exit(FatalError);
178 }
179
180 Info<< "lambda : " << lambda << nl
181 << "mu : " << mu << nl
182 << "Iters : " << iters << nl
183 << "Reading surface from " << surfFileName << " ..." << endl;
184
185 meshedSurface surf1(surfFileName);
186
187 Info<< "Faces : " << surf1.size() << nl
188 << "Vertices : " << surf1.nPoints() << nl
189 << "Bounding Box: " << boundBox(surf1.localPoints()) << endl;
190
191 bitSet fixedPoints(surf1.localPoints().size(), false);
192
193 if (args.found("featureFile"))
194 {
195 const auto featureFileName = args.get<fileName>("featureFile");
196 Info<< "Reading features from " << featureFileName << " ..." << endl;
197
198 edgeMesh feMesh(featureFileName);
199
200 getFixedPoints(feMesh, surf1.localPoints(), fixedPoints);
201
202 Info<< "Number of fixed points on surface = " << fixedPoints.count()
203 << endl;
204 }
205
206 for (label iter = 0; iter < iters; iter++)
207 {
208 // Lambda
209 {
210 pointField newLocalPoints
211 (
212 (1 - lambda)*surf1.localPoints()
213 + lambda*avg(surf1, fixedPoints)
214 );
215
216 pointField newPoints(surf1.points());
217 UIndirectList<point>(newPoints, surf1.meshPoints()) =
218 newLocalPoints;
219
220 surf1.movePoints(newPoints);
221 }
222
223 // Mu
224 if (mu != 0)
225 {
226 pointField newLocalPoints
227 (
228 (1 + mu)*surf1.localPoints()
229 - mu*avg(surf1, fixedPoints)
230 );
231
232 pointField newPoints(surf1.points());
233 UIndirectList<point>(newPoints, surf1.meshPoints()) =
234 newLocalPoints;
235
236 surf1.movePoints(newPoints);
237 }
238 }
239
240 Info<< "Writing surface to " << outFileName << " ..." << endl;
241 surf1.write(outFileName);
242
243 Info<< "End\n" << endl;
244
245 return 0;
246}
247
248
249// ************************************************************************* //
virtual bool write(const token &tok)=0
Write token to stream or otherwise handle it.
A List with indirect addressing. Like IndirectList but does not store addressing.
Definition: IndirectList.H:79
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 found(const word &optName) const
Return true if the named option is found.
Definition: argListI.H:178
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
Definition: bitSet.H:66
unsigned int count(const bool on=true) const
Count number of bits set.
Definition: bitSetI.H:500
void set(const bitSet &bitset)
Set specified bits from another bitset.
Definition: bitSetI.H:590
bool test(const label pos) const
Test value at specified position, never auto-vivify entries.
Definition: bitSetI.H:521
A bounding box defined in terms of min/max extrema points.
Definition: boundBox.H:64
Mesh data needed to do the Finite Area discretisation.
Definition: edgeFaMesh.H:56
const pointField & points() const noexcept
Return points.
Definition: edgeMeshI.H:99
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:66
A class for handling file names.
Definition: fileName.H:76
A class for managing temporary objects.
Definition: tmp.H:65
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
const volScalarField & mu
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
const pointField & points
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
Determine correspondence between points. See below.
#define WarningInFunction
Report a warning using Foam::Warning.
Namespace for OpenFOAM.
messageStream Info
Information stream (stdout output on master, null elsewhere)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
bool matchPoints(const UList< point > &pts0, const UList< point > &pts1, const UList< scalar > &matchDistance, const bool verbose, labelList &from0To1, const point &origin=point::zero)
Determine correspondence between pointFields. Gets passed.
Definition: matchPoints.C:36
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
loopControl iters(runTime, aMesh.solutionDict(), "solution")
Foam::argList args(argc, argv)
volScalarField & e
Definition: createFields.H:11
dimensionedScalar lambda("lambda", dimTime/sqr(dimLength), laminarTransport)
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333