60 if (!
mesh.isInternalFace(facei))
63 for (
const label edgei: faceEdges)
65 const labelList& edgeIFaces = edgeFaces[edgei];
66 for (
const label neiFacei : edgeIFaces)
68 if (neiFacei != facei && !
mesh.isInternalFace(neiFacei))
71 mesh.boundaryMesh().whichPatch(neiFacei);
72 if (!isA<emptyFvPatch>(
mesh.boundary()[patchi]))
87 const cutFaceIso& cutFace,
88 DynamicList<vector>& isoSurfPts,
89 DynamicList<face>& isoSurfFaces,
91 List<DynamicList<label>>& cuttingFacesPerMeshFace
94 const fvMesh&
mesh = zones_.mesh();
100 if (!isA<coupledFvPatch>(cutPatch) && !isA<emptyFvPatch>(cutPatch))
104 addCuttingFaceToIsoline
106 cutFace.subFacePoints(),
109 cuttingFacesPerMeshFace,
118 isoSurfFaces.size() - 1
131 const DynamicList<point>& facePoints,
132 const label nSerialPatches,
133 const DynamicList<label>& cellCutFaces,
134 const List<DynamicList<label>>& cuttingFacesPerMeshFace,
135 DynamicList<vector>& isoSurfPts,
136 DynamicList<face>& isoSurfFaces,
140 if (facePoints.size() > 2)
145 labelList uniquePointIDs(facePoints.size(), -1);
146 DynamicList<point> uniqueFacePoints(facePoints.size());
147 DynamicList<label> uniqueFacePointEdges(facePoints.size());
148 label addedPoints =
Zero;
151 bool foundInNei =
false;
152 for (
const label cutMeshFacei : cellCutFaces)
160 cuttingFacesPerMeshFace[cutMeshFacei],
173 uniquePointIDs[
pi] = isoSurfPts.
size() + addedPoints;
175 uniqueFacePoints.push_back(facePoints[
pi]);
179 face isoFace(uniquePointIDs);
180 isoSurfPts.append(uniqueFacePoints);
181 isoSurfFaces.append(isoFace);
194 const DynamicList<label>& cuttingFaces,
195 const DynamicList<point>& isoSurfPts,
196 const DynamicList<face>& isoSurfFaces,
200 for (
const label cuttingFacei : cuttingFaces)
202 const face& cuttingFace = isoSurfFaces[cuttingFacei];
203 for (
const label neiPi : cuttingFace)
205 if (
mag(pointi - isoSurfPts[neiPi]) < SMALL)
207 uniquePointIDs[pointID] = neiPi;
219 const Map<label>& addedFaces,
220 const scalar isoValue,
221 DynamicList<vector>& isoSurfPts,
222 DynamicList<face>& isoSurfFaces,
224 label& nChangedFaces,
226 List<wallPointData<label>>& changedFacesInfo,
228 List<DynamicList<label>>& cuttingFacesPerMeshFace
231 const fvMesh&
mesh = zones_.mesh();
237 bool isWall = isA<wallFvPatch>(
patch);
238 if (!isA<emptyFvPatch>(
patch) && !isA<coupledFvPatch>(
patch))
240 const label start =
patch.start();
243 const label gFacei = start + facei;
249 for (
const label pti :
mesh.
faces()[gFacei])
251 isFluid = isFluid && (pointY[pti] >= isoValue);
253 if (isFluid && !addedFaces.found(gFacei))
259 meshFaceToChangedFace_.insert(gFacei, nChangedFaces);
261 changedFacesInfo[nChangedFaces] =
268 changedFaces[nChangedFaces] = gFacei;
270 changedFaceToCutFace.push_back(isoSurfFaces.size());
276 addCuttingFaceToIsoline
281 cuttingFacesPerMeshFace,
288 cuttingFacesPerMeshFace[gFacei].push_back
290 isoSurfFaces.size() - 1
305 const label nSerialPatches
308 const fvMesh&
mesh = zones_.mesh();
312 fileName localName(
"topOIsoSurface" +
timeName);
313 fileName fname(isoSurfFolder_/localName);
314 vtk::surfaceWriter vtkFile(
pts, faces, fname);
316 vtkFile.writeGeometry();
317 vtkFile.beginCellData(1);
318 vtkFile.writeCellData(
"zoneIds", zoneIds);
329 autoPtr<meshedSurf> surf(
nullptr);
334 autoPtr<meshedSurf>::NewFrom<mergedSurf>
344 autoPtr<meshedSurf>::NewFrom<meshedSurfRef>
346 pts, faces, zoneIds, faceIds
353 const faceList& surfFaces = surf().faces();
354 const labelList& surfZoneIds = surf().zoneIds();
357 labelList zoneSizes(nSerialPatches + 1, 0);
360 ++zoneSizes[surfZoneIds[fi]];
364 labelList cumulZoneSizes(nSerialPatches + 1, 0);
365 for (label zi = 1; zi < cumulZoneSizes.size(); ++zi)
367 cumulZoneSizes[zi] = cumulZoneSizes[zi - 1] + zoneSizes[zi - 1];
374 labelList passedFacesPerZone(surfZoneIds.size(), 0);
378 const label zi = surfZoneIds[fi];
379 faceMap[cumulZoneSizes[zi] + passedFacesPerZone[zi]] = fi;
380 ++passedFacesPerZone[zi];
385 List<surfZone> surfZones(nSerialPatches + 1);
391 if (!isA<coupledFvPatch>(
patch) && !isA<emptyFvPatch>(
patch))
394 surfZone(
name, zoneSizes[zi], cumulZoneSizes[zi], zi);
402 zoneSizes[nSerialPatches],
403 cumulZoneSizes[nSerialPatches],
406 surfZones.setSize(zi + 1);
408 MeshedSurfaceProxy<face>
430 Foam::topOVariablesBase::topOVariablesBase
433 const dictionary&
dict 444 IOobject::READ_IF_PRESENT,
451 (
mesh.time().globalPath()/
"optimisation"/
"topOIsoSurfaces"),
452 meshFaceToChangedFace_(),
468 const word& interpolationFieldName
483 const word& designVariablesName,
484 const word& interpolationFieldName
487 if (designVariablesName ==
"beta")
499 const scalar isoValue,
501 List<wallPointData<label>>& changedFacesInfo
504 const fvMesh&
mesh = zones_.mesh();
507 label nChangedFaces(0);
513 pointY.rename(
"pointY");
518 cutCellIso cutCell(
mesh, pointY);
519 cutFaceIso cutFace(
mesh, pointY);
520 DynamicList<face> isoSurfFaces(
mesh.
nFaces()/100);
528 meshFaceToChangedFace_.clearStorage();
532 DynamicList<label> changedFaceToCutFace;
535 List<DynamicList<label>> cuttingFacesPerMeshFace(
mesh.
nFaces());
542 if (!cutCell.calcSubCell(celli, isoValue))
544 const vector& cuttingCf = cutCell.faceCentre();
545 vector cuttingNf = cutCell.faceArea();
547 DynamicList<label> cellCutFaces(
cells[celli].size());
553 const label facei =
cells[celli][fi];
557 const scalar distSqr =
558 magSqr((Cf - cuttingCf) & cuttingNf);
560 Cf + ((cuttingCf - Cf) & cuttingNf)*cuttingNf;
561 cellCutFaces.push_back(facei);
563 if (!meshFaceToChangedFace_.found(facei))
567 bool addedToIsoSurf = addCutBoundaryFaceToIsoline
574 cuttingFacesPerMeshFace
579 meshFaceToChangedFace_.insert(facei, nChangedFaces);
581 changedFacesInfo[nChangedFaces] =
588 changedFaces[nChangedFaces] = facei;
590 changedFaceToCutFace.push_back(isoSurfFaces.size());
596 label visitedFace = meshFaceToChangedFace_.at(facei);
597 if (distSqr < changedFacesInfo[visitedFace].distSqr())
600 changedFacesInfo[visitedFace] =
608 changedFaceToCutFace[visitedFace] =
617 addCuttingFaceToIsoline
619 cutCell.facePoints(),
622 cuttingFacesPerMeshFace,
629 for (
const label facei : cellCutFaces)
631 cuttingFacesPerMeshFace[facei].push_back
633 isoSurfFaces.size() - 1
642 addBoundaryFacesToIsoline
645 meshFaceToChangedFace_,
653 changedFaceToCutFace,
654 cuttingFacesPerMeshFace
657 changedFaces.setSize(nChangedFaces);
658 changedFacesInfo.setSize(nChangedFaces);
661 surfPoints_.transfer(isoSurfPts);
662 surfFaces_.transfer(isoSurfFaces);
665 writeSurfaceFiles(surfPoints_, surfFaces_, zoneIds, nSerialPatches);
673 const label changedFacesOffset =
676 for (
auto& info : changedFacesInfo)
678 info.data() += changedFacesOffset;
bool addCuttingFaceToIsoline(const DynamicList< point > &facePoints, const label nSerialPatches, const DynamicList< label > &cellCutFaces, const List< DynamicList< label >> &cuttingFacesPerMeshFace, DynamicList< vector > &isoSurfPts, DynamicList< face > &isoSurfFaces, DynamicList< label > &zoneIDs) const
Add the cutting face to the zero level-set iso-surface.
void size(const label n)
Older name for setAddressableSize.
const labelIOList & zoneIDs
List< cell > cellList
List of cell.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
label nPoints() const noexcept
Number of mesh points.
const surfaceVectorField & Cf() const
Return face centres as surfaceVectorField.
void append(const T &val)
Append an element at the end of the list.
bool interpolate(const vector &p1, const vector &p2, const vector &o, vector &n, scalar l)
topOZones zones_
Cell zones useful for defining the constant and changing parts of the domain in topO.
static bool & parRun() noexcept
Test if this a parallel run.
const cellList & cells() const
Vector< Cmpt > & normalise(const scalar tol=ROOTVSMALL)
Inplace normalise the vector by its magnitude.
bool isDuplicatePoint(const label pointID, const vector &pointi, const DynamicList< label > &cuttingFaces, const DynamicList< point > &isoSurfPts, const DynamicList< face > &isoSurfFaces, labelList &uniquePointIDs) const
bool writeTime() const noexcept
True if this is a write interval.
const Time & time() const
Return the top-level database.
label nFaces() const noexcept
Number of mesh faces.
label calcSubFace(const label faceI, const scalar cutValue)
Calculate cut points along edges of faceI.
Macros for easy insertion into run-time selection tables.
bool isInternalFace(const label faceIndex) const noexcept
Return true if given face label is internal to the mesh.
Pair< int > faceMap(const label facePi, const face &faceP, const label faceNi, const face &faceN)
DynamicList< label > faceFaces(const label facei) const
Construct facesFaces for a given boundary face.
virtual const pointField & points() const
Return raw points.
#define forAll(list, i)
Loop across all elements in list.
GeometricField< scalar, fvPatchField, volMesh > volScalarField
static scalar defaultMergeDim
The default merge dimension (1e-8)
List< face > faceList
List of faces.
void write(vtk::formatter &fmt, const Type &val, const label n=1)
Component-wise write of a value (N times)
vectorField pointField
pointField is a vectorField.
bool mkDir(const fileName &pathName, mode_t mode=0777)
Make a directory and return an error if it could not be created.
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for expressions::valueTypeCode::INVALID.
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects...
virtual void sourceTerm(DimensionedField< scalar, volMesh > &field, const topOInterpolationFunction &interpolationFunc, const scalar betaMax, const word &interpolationFieldName="beta") const
Populate source terms for the flow equations.
const polyBoundaryMesh & boundaryMesh() const noexcept
Return boundary mesh.
A class for handling words, derived from Foam::string.
fileName isoSurfFolder_
Folder name holding the zero level-set iso-surface.
virtual void sourceTermSensitivities(scalarField &sens, const topOInterpolationFunction &interpolationFunc, const scalar betaMax, const word &designVariablesName, const word &interpolationFieldName="beta") const
Post-processing sensitivities due to interpolations based on the indicator fields.
constexpr scalar pi(M_PI)
virtual const faceList & faces() const
Return raw faces.
label whichPatch(const label meshFacei) const
Return patch index for a given mesh face index. Uses binary search.
GeometricField< scalar, pointPatchField, pointMesh > pointScalarField
virtual tmp< scalarField > derivative(const scalarField &arg) const =0
Return of function with respect to the argument field.
static word timeName(const scalar t, const int precision=precision_)
Return a time name for the given scalar time value formatted with the given precision.
int debug
Static debugging option.
Base class for selecting the betaMax value, i.e. the value multiplying the Brinkman penalisation term...
defineTypeNameAndDebug(combustionModel, 0)
void writeFluidSolidInterface(const volScalarField &indicator, const scalar isoValue, labelList &changedFaces, List< wallPointData< label >> &changedFacesInfo)
Write the fluid-solid interface to files.
const fvMesh & mesh() const
Const reference to mesh.
void push_back(const T &val)
Copy append an element to the end of this list.
void addBoundaryFacesToIsoline(const pointScalarField &pointY, const Map< label > &addedFaces, const scalar isoValue, DynamicList< vector > &isoSurfPts, DynamicList< face > &isoSurfFaces, DynamicList< label > &zoneIDs, label &nChangedFaces, labelList &changedFaces, List< wallPointData< label >> &changedFacesInfo, labelList &changedFaceToCutFace, List< DynamicList< label >> &cuttingFacesPerMeshFace)
Check whether the boundary faces of the initial domain belong to the fluid part and add them to the s...
Mesh data needed to do the Finite Volume discretisation.
static bool master(const label communicator=worldComm)
True if process corresponds to the master rank in the communicator.
Field with dimensions and associated with geometry type GeoMesh which is used to size the field and a...
void push_back(T *ptr)
Append an element to the end of the list.
const std::string patch
OpenFOAM patch number as a std::string.
virtual void interpolate(const scalarField &arg, scalarField &res) const =0
Interpolate argument to result.
const fvBoundaryMesh & boundary() const noexcept
Return reference to boundary mesh.
dimensionedScalar beta("beta", dimless/dimTemperature, laminarTransport)
List< label > labelList
A List of labels.
bool addCutBoundaryFaceToIsoline(const label facei, const cutFaceIso &cutFace, DynamicList< vector > &isoSurfPts, DynamicList< face > &isoSurfFaces, DynamicList< label > &zoneIDs, List< DynamicList< label >> &cuttingFacesPerMeshFace) const
Check whether the cutFace intersects the boundary of the initial domain and add fluid part of the int...
label nNonProcessor() const
The number of patches before the first processor patch.
void writeSurfaceFiles(const pointField &pts, const faceList &faces, const labelList &zoneIds, const label nSerialPatches) const
Write the surface describing the fluid domain to stl and vtp files.
Base class for all design variables related to topology optimisation (topO). Provides the lookup func...
static label calcOffset(const label localSize, const label comm=UPstream::worldComm, const bool checkOverflow=false)
Calculate globally-consistent local start offset based on the local input size(s).
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
static constexpr const zero Zero
Global zero (0)