Commit 2aba8cff authored by Franjo's avatar Franjo

Initial implementation of the whole methodology

parent ace1d529
......@@ -27,20 +27,10 @@ Description
#include "demandDrivenData.H"
#include "symmetryPlaneOptimisation.H"
#include "partTetMesh.H"
#include "meshSurfaceEngine.H"
#include "meshSurfaceEngineModifier.H"
#include "partTetMeshSimplex.H"
#include "meshUntangler.H"
#include "volumeOptimizer.H"
#include "knuppMetric.H"
#include <map>
# ifdef USE_OMP
#include <omp.h>
# endif
#include "polyMeshGenAddressing.H"
#include "helperFunctions.H"
#include "polyMeshGenChecks.H"
#include "meshOptimizer.H"
// #define DEBUGSearch
......@@ -64,6 +54,7 @@ void symmetryPlaneOptimisation::detectSymmetryPlanes()
forAll(boundaries, patchI)
{
Info << "Patch " << patchI << " is of type " << boundaries[patchI].type() << endl;
if( boundaries[patchI].type() == "symmetryPlane" )
{
std::pair<vector, label>& cs = centreSum[patchI];
......@@ -108,14 +99,82 @@ void symmetryPlaneOptimisation::detectSymmetryPlanes()
const std::pair<vector, label>& ns = normalSum[it->first];
const point n = ns.first / ns.second;
const word pName = boundaries[it->first].patchName();
symmetryPlanes_.insert(std::make_pair(pName, plane(c, n)));
symmetryPlanes_.insert(std::make_pair(it->first, plane(c, n)));
}
}
void symmetryPlaneOptimisation::pointInPlanes(VRWGraph&) const
void symmetryPlaneOptimisation::pointInPlanes(VRWGraph& pointInPlanes) const
{
const PtrList<boundaryPatch>& boundaries = mesh_.boundaries();
const pointFieldPMG& points = mesh_.points();
const faceListPMG& faces = mesh_.faces();
pointInPlanes.clear();
pointInPlanes.setSize(points.size());
forAll(boundaries, patchI)
{
if( boundaries[patchI].type() == "symmetryPlane" )
{
const label start = boundaries[patchI].patchStart();
const label end = start + boundaries[patchI].patchSize();
for(label faceI=start;faceI<end;++faceI)
{
const face& f = faces[faceI];
forAll(f, pI)
pointInPlanes.appendIfNotIn(f[pI], patchI);
}
}
}
if( Pstream::parRun() )
{
const Map<label>& globalToLocal =
mesh_.addressingData().globalToLocalPointAddressing();
const VRWGraph& pointAtProcs = mesh_.addressingData().pointAtProcs();
const DynList<label>& neiProcs =
mesh_.addressingData().pointNeiProcs();
std::map<label, labelLongList> exchangeData;
forAll(neiProcs, i)
exchangeData[neiProcs[i]].clear();
forAllConstIter(Map<label>, globalToLocal, it)
{
const label pointI = it();
if( pointInPlanes.sizeOfRow(pointI) == 0 )
continue;
forAllRow(pointAtProcs, pointI, i)
{
const label neiProc = pointAtProcs(pointI, i);
if( neiProc == Pstream::myProcNo() )
continue;
labelLongList& dataToSend = exchangeData[neiProc];
dataToSend.append(it.key());
dataToSend.append(pointInPlanes.sizeOfRow(pointI));
forAllRow(pointInPlanes, pointI, pipI)
dataToSend.append(pointInPlanes(pointI, pipI));
}
}
labelLongList receivedData;
help::exchangeMap(exchangeData, receivedData);
for(label counter=0;counter<receivedData.size();)
{
const label pointI = globalToLocal[receivedData[counter++]];
const label size = receivedData[counter++];
for(label i=0;i<size;++i)
pointInPlanes.appendIfNotIn(pointI, receivedData[counter++]);
}
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
......@@ -137,7 +196,55 @@ symmetryPlaneOptimisation::~symmetryPlaneOptimisation()
void symmetryPlaneOptimisation::optimizeSymmetryPlanes()
{
pointFieldPMG& points = mesh_.points();
Info << "Detected " << label(symmetryPlanes_.size()) << " symmetry planes in the mesh" << endl;
VRWGraph pointInPlane;
pointInPlanes(pointInPlane);
forAll(pointInPlane, pointI)
{
const label nPlanes = pointInPlane.sizeOfRow(pointI);
if( nPlanes > 3 )
{
WarningIn
(
"void symmetryPlaneOptimisation::optimizeSymmetryPlanes()"
) << "Point " << pointI << " is in more than three symmetry"
<< " planes. Cannot move it" << endl;
}
point& p = points[pointI];
vector disp(vector::zero);
for(label plI=0;plI<nPlanes;++plI)
{
//- point is in a plane
std::map<label, plane>::const_iterator it =
symmetryPlanes_.find(pointInPlane(pointI, 0));
const point newP = it->second.nearestPoint(points[pointI]);
disp += newP - p;
}
p += disp;
}
labelHashSet badFaces;
polyMeshGenChecks::checkFacePyramids(mesh_, false, VSMALL, &badFaces);
if( badFaces.size() )
{
WarningIn
(
"void symmetryPlaneOptimisation::optimizeSymmetryPlanes()"
) << "Bad quality or inverted faces found in the mesh" << endl;
const label badFacesId = mesh_.addFaceSubset("invalidFaces");
forAllConstIter(labelHashSet, badFaces, it)
mesh_.addFaceToSubset(badFacesId, it.key());
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......
......@@ -61,7 +61,7 @@ class symmetryPlaneOptimisation
polyMeshGen& mesh_;
//- symmetry planes in the mesh
std::map<word, plane> symmetryPlanes_;
std::map<label, plane> symmetryPlanes_;
// Private member functions
//- detect symmetry planes
......
......@@ -44,8 +44,10 @@ int main(int argc, char *argv[])
polyMeshGen pmg(runTime);
pmg.read();
Info << "Starting optimisation of symmetry planes" << endl;
symmetryPlaneOptimisation(pmg).optimizeSymmetryPlanes();
Info << "Writing mesh" << endl;
pmg.write();
Info << "End\n" << endl;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment