controlMeshRefinement.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) 2013-2015 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
30#include "OFstream.H"
31
32// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
33
34namespace Foam
35{
36 defineTypeNameAndDebug(controlMeshRefinement, 0);
37}
38
39// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
40
41
42// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
43
44Foam::scalar Foam::controlMeshRefinement::calcFirstDerivative
45(
46 const Foam::point& a,
47 const scalar& cellSizeA,
48 const Foam::point& b,
49 const scalar& cellSizeB
50) const
51{
52 return (cellSizeA - cellSizeB)/mag(a - b);
53}
54
55
56//Foam::scalar Foam::controlMeshRefinement::calcSecondDerivative
57//(
58// const Foam::point& a,
59// const scalar& cellSizeA,
60// const Foam::point& midPoint,
61// const scalar& cellSizeMid,
62// const Foam::point& b,
63// const scalar& cellSizeB
64//) const
65//{
66// return (cellSizeA - 2*cellSizeMid + cellSizeB)/magSqr((a - b)/2);
67//}
68
69
70bool Foam::controlMeshRefinement::detectEdge
71(
72 const Foam::point& startPt,
73 const Foam::point& endPt,
74 pointHit& pointFound,
75 const scalar tolSqr,
76 const scalar secondDerivTolSqr
77) const
78{
79 Foam::point a(startPt);
80 Foam::point b(endPt);
81
82 Foam::point midPoint = (a + b)/2.0;
83
84 label nIterations = 0;
85
86 while (true)
87 {
88 nIterations++;
89
90 if
91 (
92 magSqr(a - b) < tolSqr
93 )
94 {
95 pointFound.setPoint(midPoint);
96 pointFound.setHit();
97
98 return true;
99 }
100
101 // Split into two regions
102
103 scalar cellSizeA = sizeControls_.cellSize(a);
104 scalar cellSizeB = sizeControls_.cellSize(b);
105
106// if (magSqr(cellSizeA - cellSizeB) < 1e-6)
107// {
108// return false;
109// }
110
111 scalar cellSizeMid = sizeControls_.cellSize(midPoint);
112
113 // Region 1
114 Foam::point midPoint1 = (a + midPoint)/2.0;
115 const scalar cellSizeMid1 = sizeControls_.cellSize(midPoint1);
116
117// scalar firstDerivative1 =
118// calcFirstDerivative(cellSizeA, cellSizeMid);
119
120 scalar secondDerivative1 =
121 calcSecondDerivative
122 (
123 a,
124 cellSizeA,
125 midPoint1,
126 cellSizeMid1,
127 midPoint,
128 cellSizeMid
129 );
130
131 // Region 2
132 Foam::point midPoint2 = (midPoint + b)/2.0;
133 const scalar cellSizeMid2 = sizeControls_.cellSize(midPoint2);
134
135// scalar firstDerivative2 =
136// calcFirstDerivative(f, cellSizeMid, cellSizeB);
137
138 scalar secondDerivative2 =
139 calcSecondDerivative
140 (
141 midPoint,
142 cellSizeMid,
143 midPoint2,
144 cellSizeMid2,
145 b,
146 cellSizeB
147 );
148
149 // Neither region appears to have an inflection
150 // To be sure should use higher order derivatives
151 if
152 (
153 magSqr(secondDerivative1) < secondDerivTolSqr
154 && magSqr(secondDerivative2) < secondDerivTolSqr
155 )
156 {
157 return false;
158 }
159
160 // Pick region with greatest second derivative
161 if (magSqr(secondDerivative1) > magSqr(secondDerivative2))
162 {
163 b = midPoint;
164 midPoint = midPoint1;
165 }
166 else
167 {
168 a = midPoint;
169 midPoint = midPoint2;
170 }
171 }
172}
173
174
175Foam::pointHit Foam::controlMeshRefinement::findDiscontinuities
176(
177 const linePointRef& l
178) const
179{
181
182 const scalar tolSqr = sqr(1e-3);
183 const scalar secondDerivTolSqr = sqr(1e-3);
184
185 detectEdge
186 (
187 l.start(),
188 l.end(),
189 p,
190 tolSqr,
191 secondDerivTolSqr
192 );
193
194 return p;
195}
196
197
198// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
199
200
201// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
202
204(
205 cellShapeControl& shapeController
206)
207:
208 shapeController_(shapeController),
209 mesh_(shapeController.shapeControlMesh()),
210 sizeControls_(shapeController.sizeAndAlignment()),
211 geometryToConformTo_(sizeControls_.geometryToConformTo())
212{}
213
214
215// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
216
218{}
219
220
221// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
222
224(
225 const autoPtr<backgroundMeshDecomposition>& decomposition
226)
227{
228 if (shapeController_.shapeControlMesh().vertexCount() > 0)
229 {
230 // Mesh already populated.
231 Info<< "Cell size and alignment mesh already populated." << endl;
232 return;
233 }
234
235 autoPtr<boundBox> overallBoundBox;
236
237 // Need to pass in the background mesh decomposition so that can test if
238 // a point to insert is on the processor.
239 if (Pstream::parRun())
240 {
241// overallBoundBox.set(new boundBox(decomposition().procBounds()));
242 }
243 else
244 {
245// overallBoundBox.set
246// (
247// new boundBox(geometryToConformTo_.geometry().bounds())
248// );
249//
250// mesh_.insertBoundingPoints
251// (
252// overallBoundBox(),
253// sizeControls_
254// );
255 }
256
257 Map<label> priorityMap;
258
259 const PtrList<cellSizeAndAlignmentControl>& controlFunctions =
260 sizeControls_.controlFunctions();
261
262 forAll(controlFunctions, fI)
263 {
264 const cellSizeAndAlignmentControl& controlFunction =
265 controlFunctions[fI];
266
267 const Switch& forceInsertion =
268 controlFunction.forceInitialPointInsertion();
269
270 Info<< "Inserting points from " << controlFunction.name()
271 << " (" << controlFunction.type() << ")" << endl;
272 Info<< " Force insertion is " << forceInsertion.c_str() << endl;
273
274 pointField pts;
275 scalarField sizes;
276 triadField alignments;
277
278 controlFunction.initialVertices(pts, sizes, alignments);
279
280 Info<< " Got initial vertices list of size " << pts.size() << endl;
281
282 List<Vb> vertices(pts.size());
283
284 // Clip the minimum size
285 for (label vI = 0; vI < pts.size(); ++vI)
286 {
287 vertices[vI] = Vb(pts[vI], Vb::vtInternalNearBoundary);
288
289 label maxPriority = -1;
290 scalar size = sizeControls_.cellSize(pts[vI], maxPriority);
291
292 if (maxPriority > controlFunction.maxPriority())
293 {
294 vertices[vI].targetCellSize() = max
295 (
296 size,
297 shapeController_.minimumCellSize()
298 );
299 }
300// else if (maxPriority == controlFunction.maxPriority())
301// {
302// vertices[vI].targetCellSize() = max
303// (
304// min(sizes[vI], size),
305// shapeController_.minimumCellSize()
306// );
307// }
308 else
309 {
310 vertices[vI].targetCellSize() = max
311 (
312 sizes[vI],
313 shapeController_.minimumCellSize()
314 );
315 }
316
317 vertices[vI].alignment() = alignments[vI];
318 }
319
320 Info<< " Clipped minimum size" << endl;
321
322 pts.clear();
323 sizes.clear();
324 alignments.clear();
325
326 bitSet keepVertex(vertices.size(), true);
327
328 forAll(vertices, vI)
329 {
330 bool keep = true;
331
333
334 if (Pstream::parRun())
335 {
336 keep = decomposition().positionOnThisProcessor(pt);
337 }
338
339 if (keep && geometryToConformTo_.wellOutside(pt, SMALL))
340 {
341 keep = false;
342 }
343
344 if (!keep)
345 {
346 keepVertex.unset(vI);
347 }
348 }
349
350 inplaceSubset(keepVertex, vertices);
351
352 const label preInsertedSize = mesh_.number_of_vertices();
353
354 Info<< " Check sizes" << endl;
355
356 forAll(vertices, vI)
357 {
358 bool insertPoint = false;
359
361
362 if
363 (
364 mesh_.dimension() < 3
365 || mesh_.is_infinite
366 (
367 mesh_.locate(vertices[vI].point())
368 )
369 )
370 {
371 insertPoint = true;
372 }
373
374 const scalar interpolatedCellSize = shapeController_.cellSize(pt);
375 const triad interpolatedAlignment =
376 shapeController_.cellAlignment(pt);
377 const scalar calculatedCellSize = vertices[vI].targetCellSize();
378 const triad calculatedAlignment = vertices[vI].alignment();
379
380 if (debug)
381 {
382 Info<< "Point = " << pt << nl
383 << " Size(interp) = " << interpolatedCellSize << nl
384 << " Size(calc) = " << calculatedCellSize << nl
385 << " Align(interp) = " << interpolatedAlignment << nl
386 << " Align(calc) = " << calculatedAlignment << nl
387 << endl;
388 }
389
390 const scalar sizeDiff =
391 mag(interpolatedCellSize - calculatedCellSize);
392 const scalar alignmentDiff =
393 diff(interpolatedAlignment, calculatedAlignment);
394
395 if (debug)
396 {
397 Info<< " size difference = " << sizeDiff << nl
398 << ", alignment difference = " << alignmentDiff << endl;
399 }
400
401 // TODO: Also need to base it on the alignments
402 if
403 (
404 sizeDiff/interpolatedCellSize > 0.1
405 || alignmentDiff > 0.15
406 )
407 {
408 insertPoint = true;
409 }
410
411 if (forceInsertion || insertPoint)
412 {
413 const label oldSize = mesh_.vertexCount();
414
415 cellShapeControlMesh::Vertex_handle insertedVert = mesh_.insert
416 (
417 pt,
418 calculatedCellSize,
419 vertices[vI].alignment(),
421 );
422
423 if (oldSize == mesh_.vertexCount() - 1)
424 {
425 priorityMap.insert
426 (
427 insertedVert->index(),
428 controlFunction.maxPriority()
429 );
430 }
431 }
432 }
433
434 //mesh_.rangeInsertWithInfo(vertices.begin(), vertices.end());
435
436 Info<< " Inserted "
437 << returnReduce
438 (
439 label(mesh_.number_of_vertices()) - preInsertedSize,
440 sumOp<label>()
441 )
442 << "/" << returnReduce(vertices.size(), sumOp<label>())
443 << endl;
444 }
445
446
447
448 forAll(controlFunctions, fI)
449 {
450 const cellSizeAndAlignmentControl& controlFunction =
451 controlFunctions[fI];
452
453 const Switch& forceInsertion =
454 controlFunction.forceInitialPointInsertion();
455
456 Info<< "Inserting points from " << controlFunction.name()
457 << " (" << controlFunction.type() << ")" << endl;
458 Info<< " Force insertion is " << forceInsertion.c_str() << endl;
459
460 DynamicList<Foam::point> extraPts;
461 DynamicList<scalar> extraSizes;
462
463 controlFunction.cellSizeFunctionVertices(extraPts, extraSizes);
464
465 List<Vb> vertices(extraPts.size());
466
467 // Clip the minimum size
468 for (label vI = 0; vI < extraPts.size(); ++vI)
469 {
470 vertices[vI] = Vb(extraPts[vI], Vb::vtUnassigned);
471
472 label maxPriority = -1;
473 scalar size = sizeControls_.cellSize(extraPts[vI], maxPriority);
474
475 if (maxPriority > controlFunction.maxPriority())
476 {
477 vertices[vI].targetCellSize() = max
478 (
479 size,
480 shapeController_.minimumCellSize()
481 );
482 }
483 else if (maxPriority == controlFunction.maxPriority())
484 {
485 vertices[vI].targetCellSize() = max
486 (
487 min(extraSizes[vI], size),
488 shapeController_.minimumCellSize()
489 );
490 }
491 else
492 {
493 vertices[vI].targetCellSize() = max
494 (
495 extraSizes[vI],
496 shapeController_.minimumCellSize()
497 );
498 }
499 }
500
501 bitSet keepVertex(vertices.size(), true);
502
503 forAll(vertices, vI)
504 {
505 bool keep = true;
506
508
509 if (Pstream::parRun())
510 {
511 keep = decomposition().positionOnThisProcessor(pt);
512 }
513
514 if (keep && geometryToConformTo_.wellOutside(pt, SMALL))
515 {
516 keep = false;
517 }
518
519 if (!keep)
520 {
521 keepVertex.unset(vI);
522 }
523 }
524
525 inplaceSubset(keepVertex, vertices);
526
527 const label preInsertedSize = mesh_.number_of_vertices();
528
529 forAll(vertices, vI)
530 {
531 bool insertPoint = false;
532
534
535 if
536 (
537 mesh_.dimension() < 3
538 || mesh_.is_infinite
539 (
540 mesh_.locate(vertices[vI].point())
541 )
542 )
543 {
544 insertPoint = true;
545 }
546
547 const scalar interpolatedCellSize = shapeController_.cellSize(pt);
548 const scalar calculatedCellSize = vertices[vI].targetCellSize();
549
550 if (debug)
551 {
552 Info<< "Point = " << pt << nl
553 << " Size(interp) = " << interpolatedCellSize << nl
554 << " Size(calc) = " << calculatedCellSize << nl
555 << endl;
556 }
557
558 const scalar sizeDiff =
559 mag(interpolatedCellSize - calculatedCellSize);
560
561 if (debug)
562 {
563 Info<< " size difference = " << sizeDiff << endl;
564 }
565
566 // TODO: Also need to base it on the alignments
567 if (sizeDiff/interpolatedCellSize > 0.1)
568 {
569 insertPoint = true;
570 }
571
572 if (forceInsertion || insertPoint)
573 {
574 // Check the priority
575
576// cellShapeControlMesh::Cell_handle ch =
577// mesh_.locate(toPoint<cellShapeControlMesh::Point>(pt));
578
579// if (mesh_.is_infinite(ch))
580// {
581// continue;
582// }
583
584// const label newPtPriority = controlFunction.maxPriority();
585
586// label highestPriority = -1;
587// for (label cI = 0; cI < 4; ++cI)
588// {
589// if (mesh_.is_infinite(ch->vertex(cI)))
590// {
591// continue;
592// }
593
594// const label vertPriority =
595// priorityMap[ch->vertex(cI)->index()];
596
597// if (vertPriority > highestPriority)
598// {
599// highestPriority = vertPriority;
600// }
601// }
602
603// if (newPtPriority >= highestPriority)
604// {
605// const label oldSize = mesh_.vertexCount();
606//
607// cellShapeControlMesh::Vertex_handle insertedVert =
608 mesh_.insert
609 (
610 pt,
611 calculatedCellSize,
612 vertices[vI].alignment(),
614 );
615
616// if (oldSize == mesh_.vertexCount() - 1)
617// {
618// priorityMap.insert
619// (
620// insertedVert->index(),
621// newPtPriority
622// );
623// }
624// }
625 }
626 }
627
628 //mesh_.rangeInsertWithInfo(vertices.begin(), vertices.end());
629
630 Info<< " Inserted extra points "
631 << returnReduce
632 (
633 label(mesh_.number_of_vertices()) - preInsertedSize,
634 sumOp<label>()
635 )
636 << "/" << returnReduce(vertices.size(), sumOp<label>())
637 << endl;
638 }
639
640 // Change cell size function of bounding points to be consistent
641 // with their nearest neighbours
642// for
643// (
644// CellSizeDelaunay::Finite_vertices_iterator vit =
645// mesh_.finite_vertices_begin();
646// vit != mesh_.finite_vertices_end();
647// ++vit
648// )
649// {
650// if (vit->uninitialised())
651// {
652// // Get its adjacent vertices
653// std::list<CellSizeDelaunay::Vertex_handle> adjacentVertices;
654//
655// mesh_.adjacent_vertices
656// (
657// vit,
658// std::back_inserter(adjacentVertices)
659// );
660//
661// scalar totalCellSize = 0;
662// label nVerts = 0;
663//
664// for
665// (
666// std::list<CellSizeDelaunay::Vertex_handle>::iterator avit =
667// adjacentVertices.begin();
668// avit != adjacentVertices.end();
669// ++avit
670// )
671// {
672// if (!(*avit)->uninitialised())
673// {
674// totalCellSize += (*avit)->targetCellSize();
675// nVerts++;
676// }
677// }
678//
679// Pout<< "Changing " << vit->info();
680//
681// vit->targetCellSize() = totalCellSize/nVerts;
682// vit->type() = Vb::vtInternalNearBoundary;
683//
684// Pout<< "to " << vit->info() << endl;
685// }
686// }
687}
688
689
691(
692 const autoPtr<backgroundMeshDecomposition>& decomposition
693)
694{
695 Info<< "Iterate over "
696 << returnReduce(label(mesh_.number_of_finite_edges()), sumOp<label>())
697 << " cell size mesh edges" << endl;
698
699 DynamicList<Vb> verts(mesh_.number_of_vertices());
700
701 label count = 0;
702
703 for
704 (
705 CellSizeDelaunay::Finite_edges_iterator eit =
706 mesh_.finite_edges_begin();
707 eit != mesh_.finite_edges_end();
708 ++eit
709 )
710 {
711 if (count % 10000 == 0)
712 {
713 Info<< count << " edges, inserted " << verts.size()
714 << " Time = " << mesh_.time().elapsedCpuTime()
715 << endl;
716 }
717 count++;
718
719 CellSizeDelaunay::Cell_handle c = eit->first;
720 CellSizeDelaunay::Vertex_handle vA = c->vertex(eit->second);
721 CellSizeDelaunay::Vertex_handle vB = c->vertex(eit->third);
722
723 if
724 (
725 mesh_.is_infinite(vA)
726 || mesh_.is_infinite(vB)
727 || (vA->referred() && vB->referred())
728 || (vA->referred() && (vA->procIndex() > vB->procIndex()))
729 || (vB->referred() && (vB->procIndex() > vA->procIndex()))
730 )
731 {
732 continue;
733 }
734
735 pointFromPoint ptA(topoint(vA->point()));
736 pointFromPoint ptB(topoint(vB->point()));
737
738 linePointRef l(ptA, ptB);
739
740 const pointHit hitPt = findDiscontinuities(l);
741
742 if (hitPt.hit())
743 {
744 const Foam::point& pt = hitPt.hitPoint();
745
746 if (!geometryToConformTo_.inside(pt))
747 {
748 continue;
749 }
750
751 if (Pstream::parRun())
752 {
753 if (!decomposition().positionOnThisProcessor(pt))
754 {
755 continue;
756 }
757 }
758
759 verts.append
760 (
761 Vb
762 (
763 toPoint(pt),
765 )
766 );
767
768 verts.last().targetCellSize() = sizeControls_.cellSize(pt);
769 verts.last().alignment() = triad::unset;
770 }
771 }
772
773 mesh_.insertPoints(verts, false);
774
775 return verts.size();
776}
777
778
779// * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * * //
780
781
782// * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * * //
783
784
785// * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * * //
786
787
788// ************************************************************************* //
CGAL::indexedVertex< K > Vb
An indexed form of CGAL::Triangulation_vertex_base_3<K> used to keep track of the Delaunay vertices i...
Definition: indexedVertex.H:88
Describes the interaction of a face and a point. It carries the info of a successful hit and (if succ...
Definition: PointHit.H:54
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
static bool & parRun() noexcept
Test if this a parallel run.
Definition: UPstream.H:433
static constexpr direction size() noexcept
The number of elements in the VectorSpace = Ncmpts.
Definition: VectorSpace.H:176
CellSizeDelaunay::Vertex_handle Vertex_handle
void initialMeshPopulation(const autoPtr< backgroundMeshDecomposition > &decomposition)
~controlMeshRefinement()
Destructor.
label refineMesh(const autoPtr< backgroundMeshDecomposition > &decomposition)
static const complex max
complex (VGREAT,VGREAT)
Definition: complex.H:293
static const triad unset
Definition: triad.H:97
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
Definition: className.H:121
volScalarField & p
unsigned int count(const UList< bool > &bools, const bool val=true)
Count number of 'true' entries.
Definition: BitOps.H:78
const dimensionedScalar c
Speed of light in a vacuum.
Namespace for OpenFOAM.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
void inplaceSubset(const BoolListType &select, ListType &input, const bool invert=false)
Inplace extract elements of the input list when select is true.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
PointFrompoint toPoint(const Foam::point &p)
line< point, const point & > linePointRef
A line using referred points.
Definition: linePointRef.H:47
pointFromPoint topoint(const Point &P)
messageStream Info
Information stream (stdout output on master, null elsewhere)
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:44
vector point
Point is a vector.
Definition: point.H:43
Field< triad > triadField
Specialisation of Field<T> for triad.
Definition: triadFieldFwd.H:45
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
scalar diff(const triad &A, const triad &B)
Return a quantity of the difference between two triads.
Definition: triad.C:378
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
pointField vertices(const blockVertexList &bvl)
PointHit< point > pointHit
A PointIndexHit for 3D points.
Definition: pointHit.H:44
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce (copy) and return value.
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
volScalarField & b
Definition: createFields.H:27
volScalarField & e
Definition: createFields.H:11
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333