meshToMeshMethod.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2013-2019 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "meshToMeshMethod.H"
27 #include "tetOverlapVolume.H"
28 #include "OFstream.H"
29 #include "Time.H"
30 #include "treeBoundBox.H"
31 
32 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
33 
34 namespace Foam
35 {
36  defineTypeNameAndDebug(meshToMeshMethod, 0);
37  defineRunTimeSelectionTable(meshToMeshMethod, components);
38 }
39 
40 Foam::scalar Foam::meshToMeshMethod::tolerance_ = 1e-6;
41 
42 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
43 
45 {
46  boundBox intersectBb
47  (
48  max(src_.bounds().min(), tgt_.bounds().min()),
49  min(src_.bounds().max(), tgt_.bounds().max())
50  );
51 
52  intersectBb.inflate(0.01);
53 
54  const cellList& srcCells = src_.cells();
55  const faceList& srcFaces = src_.faces();
56  const pointField& srcPts = src_.points();
57 
58  DynamicList<label> cells(src_.nCells());
59  forAll(srcCells, srcI)
60  {
61  boundBox cellBb(srcCells[srcI].points(srcFaces, srcPts), false);
62  if (intersectBb.overlaps(cellBb))
63  {
64  cells.append(srcI);
65  }
66  }
67 
68  if (debug)
69  {
70  Pout<< "participating source mesh cells: " << cells.size() << endl;
71  }
72 
73  return move(cells);
74 }
75 
76 
78 (
79  const label srcCelli,
80  const label tgtCelli
81 ) const
82 {
83  scalar threshold = tolerance_*src_.cellVolumes()[srcCelli];
84 
85  tetOverlapVolume overlapEngine;
86 
87  treeBoundBox bbTgtCell(tgt_.points(), tgt_.cellPoints()[tgtCelli]);
88 
89  return overlapEngine.cellCellOverlapMinDecomp
90  (
91  src_,
92  srcCelli,
93  tgt_,
94  tgtCelli,
95  bbTgtCell,
96  threshold
97  );
98 }
99 
100 
102 (
103  const label srcCelli,
104  const label tgtCelli
105 ) const
106 {
107  tetOverlapVolume overlapEngine;
108 
109  treeBoundBox bbTgtCell(tgt_.points(), tgt_.cellPoints()[tgtCelli]);
110 
111  scalar vol = overlapEngine.cellCellOverlapVolumeMinDecomp
112  (
113  src_,
114  srcCelli,
115  tgt_,
116  tgtCelli,
117  bbTgtCell
118  );
119 
120  return vol;
121 }
122 
123 
125 (
126  const label celli,
127  const polyMesh& mesh,
128  const DynamicList<label>& visitedCells,
129  DynamicList<label>& nbrCellIDs
130 ) const
131 {
132  const labelList& nbrCells = mesh.cellCells()[celli];
133 
134  // filter out cells already visited from cell neighbours
135  forAll(nbrCells, i)
136  {
137  label nbrCelli = nbrCells[i];
138 
139  if
140  (
141  (findIndex(visitedCells, nbrCelli) == -1)
142  && (findIndex(nbrCellIDs, nbrCelli) == -1)
143  )
144  {
145  nbrCellIDs.append(nbrCelli);
146  }
147  }
148 }
149 
150 
152 (
153  labelListList& srcToTgtAddr,
154  scalarListList& srcToTgtWght,
155  labelListList& tgtToSrcAddr,
156  scalarListList& tgtToSrcWght
157 ) const
158 {
159  srcToTgtAddr.setSize(src_.nCells());
160  srcToTgtWght.setSize(src_.nCells());
161  tgtToSrcAddr.setSize(tgt_.nCells());
162  tgtToSrcWght.setSize(tgt_.nCells());
163 
164  if (!src_.nCells())
165  {
166  return false;
167  }
168  else if (!tgt_.nCells())
169  {
170  if (debug)
171  {
172  Pout<< "mesh interpolation: hhave " << src_.nCells() << " source "
173  << " cells but no target cells" << endl;
174  }
175 
176  return false;
177  }
178 
179  return true;
180 }
181 
182 
183 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
184 
186 (
187  const polyMesh& src,
188  const polyMesh& tgt
189 )
190 :
191  src_(src),
192  tgt_(tgt),
193  V_(0.0)
194 {
195  if (!src_.nCells() || !tgt_.nCells())
196  {
197  if (debug)
198  {
199  Pout<< "mesh interpolation: cells not on processor: Source cells = "
200  << src_.nCells() << ", target cells = " << tgt_.nCells()
201  << endl;
202  }
203  }
204 }
205 
206 
207 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
208 
210 {}
211 
212 
213 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
214 
216 (
217  const polyMesh& mesh1,
218  const polyMesh& mesh2,
219  const labelListList& mesh1ToMesh2Addr
220 ) const
221 {
222  Pout<< "Source size = " << mesh1.nCells() << endl;
223  Pout<< "Target size = " << mesh2.nCells() << endl;
224 
225  word fName("addressing_" + mesh1.name() + "_to_" + mesh2.name());
226 
227  if (Pstream::parRun())
228  {
229  fName = fName + "_proc" + Foam::name(Pstream::myProcNo());
230  }
231 
232  OFstream os(src_.time().path()/fName + ".obj");
233 
234  label vertI = 0;
235  forAll(mesh1ToMesh2Addr, i)
236  {
237  const labelList& addr = mesh1ToMesh2Addr[i];
238  forAll(addr, j)
239  {
240  label celli = addr[j];
241  const vector& c0 = mesh1.cellCentres()[i];
242 
243  const cell& c = mesh2.cells()[celli];
244  const pointField pts(c.points(mesh2.faces(), mesh2.points()));
245  forAll(pts, j)
246  {
247  const point& p = pts[j];
248  os << "v " << p.x() << ' ' << p.y() << ' ' << p.z() << nl;
249  vertI++;
250  os << "v " << c0.x() << ' ' << c0.y() << ' ' << c0.z()
251  << nl;
252  vertI++;
253  os << "l " << vertI - 1 << ' ' << vertI << nl;
254  }
255  }
256  }
257 }
258 
259 
260 // ************************************************************************* //
bool overlaps(const boundBox &) const
Overlaps/touches boundingBox?
Definition: boundBoxI.H:120
meshToMeshMethod(const polyMesh &src, const polyMesh &tgt)
Construct from source and target meshes.
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
const word & name() const
Return name.
Definition: IOobject.H:303
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
Output to file stream.
Definition: OFstream.H:82
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
Definition: UPstream.H:429
label nCells() const
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
virtual ~meshToMeshMethod()
Destructor.
const Cmpt & z() const
Definition: VectorI.H:87
A bounding box defined in terms of the points at its extremities.
Definition: boundBox.H:58
const cellList & cells() const
const dimensionedScalar & c
Speed of light in a vacuum.
void writeConnectivity(const polyMesh &mesh1, const polyMesh &mesh2, const labelListList &mesh1ToMesh2Addr) const
Write the connectivity (debugging)
const Cmpt & y() const
Definition: VectorI.H:81
virtual scalar interVol(const label srcCelli, const label tgtCelli) const
Return the intersection volume between two cells.
static scalar tolerance_
Tolerance used in volume overlap calculations.
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1131
const cellShapeList & cells
const pointField & points
A class for handling words, derived from string.
Definition: word.H:59
Calculates the overlap volume of two cells using tetrahedral decomposition.
void append(const T &)
Append an element at the end of the list.
Definition: ListI.H:178
labelList maskCells() const
Return src cell IDs for the overlap region.
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:296
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1156
const vectorField & cellCentres() const
const Cmpt & x() const
Definition: VectorI.H:75
void inflate(const scalar s)
Inflate box by factor*mag(span) in all dimensions.
Definition: boundBox.C:210
virtual bool initialise(labelListList &srcToTgtAddr, scalarListList &srcToTgtWght, labelListList &tgtToTgtAddr, scalarListList &tgtToTgtWght) const
static const char nl
Definition: Ostream.H:260
defineRunTimeSelectionTable(reactionRateFlameArea, dictionary)
defineTypeNameAndDebug(combustionModel, 0)
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurrence of given element and return index,.
virtual void appendNbrCells(const label tgtCelli, const polyMesh &mesh, const DynamicList< label > &visitedTgtCells, DynamicList< label > &nbrTgtCellIDs) const
Append target cell neighbour cells to cellIDs list.
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
void setSize(const label)
Reset size of List.
Definition: List.C:281
static bool & parRun()
Is this a parallel run?
Definition: UPstream.H:399
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:56
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
pointField points(const faceUList &, const pointField &) const
Return the cell vertices.
Definition: cell.C:102
Standard boundBox + extra functionality for use in octree.
Definition: treeBoundBox.H:87
bool cellCellOverlapMinDecomp(const primitiveMesh &meshA, const label cellAI, const primitiveMesh &meshB, const label cellBI, const treeBoundBox &cellBbB, const scalar threshold=0.0) const
Return true if olverlap volume is greater than threshold.
const doubleScalar e
Elementary charge.
Definition: doubleScalar.H:105
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
volScalarField & p
virtual bool intersect(const label srcCelli, const label tgtCelli) const
Return the true if cells intersect.
const labelListList & cellCells() const
scalar cellCellOverlapVolumeMinDecomp(const primitiveMesh &meshA, const label cellAI, const primitiveMesh &meshB, const label cellBI, const treeBoundBox &cellBbB) const
Calculates the overlap volume.
Namespace for OpenFOAM.