From 86f9541c12a8e751980e647475130fe707bdd837 Mon Sep 17 00:00:00 2001 From: "Oliver Borm (boroli)" Date: Tue, 8 Mar 2011 17:42:21 +0100 Subject: [PATCH 1/8] mesquite hex and prism support --- .../mesquiteMotionSolver.C | 343 +++++++++++++++--- .../mesquiteMotionSolver.H | 5 +- 2 files changed, 295 insertions(+), 53 deletions(-) diff --git a/src/dynamicMesh/meshMotion/mesquiteMotionSolver/mesquiteMotionSolver.C b/src/dynamicMesh/meshMotion/mesquiteMotionSolver/mesquiteMotionSolver.C index f8a13b105..40b3c0218 100644 --- a/src/dynamicMesh/meshMotion/mesquiteMotionSolver/mesquiteMotionSolver.C +++ b/src/dynamicMesh/meshMotion/mesquiteMotionSolver/mesquiteMotionSolver.C @@ -36,6 +36,13 @@ License #include "pointPatchField.H" #include "pointMesh.H" +#include "hexMatcher.H" +#include "tetMatcher.H" +#include "prismMatcher.H" +#include "pyrMatcher.H" +// #include "wedgeMatcher.H" +// #include "tetWedgeMatcher.H" + // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // namespace Foam @@ -73,6 +80,7 @@ mesquiteMotionSolver::mesquiteMotionSolver vtxCoords_(NULL), cellToNode_(NULL), fixFlags_(NULL), + mixedTypes_(NULL), refPoints_ ( IOobject @@ -118,6 +126,7 @@ mesquiteMotionSolver::mesquiteMotionSolver vtxCoords_(NULL), cellToNode_(NULL), fixFlags_(NULL), + mixedTypes_(NULL), refPoints_ ( IOobject @@ -153,11 +162,13 @@ void mesquiteMotionSolver::clearOut() delete [] vtxCoords_; delete [] cellToNode_; delete [] fixFlags_; + delete [] mixedTypes_; // Reset to NULL vtxCoords_ = NULL; cellToNode_ = NULL; fixFlags_ = NULL; + mixedTypes_ = NULL; } // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // @@ -360,13 +371,87 @@ void mesquiteMotionSolver::readOptions() ) ); - // Read the optimization metric - qMetric_ = word(optionsDict.lookup("optMetric")); + qMetricTable_.insert + ( + "EdgeLengthQuality", + autoPtr + ( + new Mesquite::EdgeLengthQualityMetric + ) + ); - if (!qMetricTable_.found(qMetric_)) + qMetricTable_.insert + ( + "EdgeLength", + autoPtr + ( + new Mesquite::EdgeLengthMetric + ) + ); + + qMetricTable_.insert + ( + "UntangleBeta", + autoPtr + ( + new Mesquite::UntangleBetaQualityMetric + ) + ); + + qMetricTable_.insert + ( + "LocalSize", + autoPtr + ( + new Mesquite::LocalSizeQualityMetric + ) + ); + +// qMetricTable_.insert +// ( +// "VolumeRatio", +// autoPtr +// ( +// new Mesquite::VolumeRatioQualityMetric +// ) +// ); + +// qMetricTable_.insert +// ( +// "Multiply", +// autoPtr +// ( +// new Mesquite::MultiplyQualityMetric +// ) +// ); + + qMetric_.setSize(1); + + if (optionsDict.found("optMetrics")) + { + const dictionary& optMetricsDict = optionsDict.subDict("optMetrics"); + + if (optMetricsDict.found("secondMetric")) + { + qMetric_.setSize(2); + qMetric_[1] = word(optMetricsDict.lookup("secondMetric")); + } + + if (optMetricsDict.found("firstMetric")) + { + qMetric_[0] = word(optMetricsDict.lookup("firstMetric")); + } + } + else + { + // Read the optimization metric + qMetric_[0] = word(optionsDict.lookup("optMetric")); + } + + if (!qMetricTable_.found(qMetric_[0])) { FatalErrorIn("void mesquiteMotionSolver::readOptions()") - << "Unrecognized quality metric: " << qMetric_ << nl + << "Unrecognized quality metric: " << qMetric_[0] << nl << "Available metrics are: " << nl << qMetricTable_.toc() << abort(FatalError); } @@ -469,7 +554,7 @@ void mesquiteMotionSolver::readOptions() ( new Mesquite::LInfTemplate ( - &qMetricTable_[qMetric_]() + &qMetricTable_[qMetric_[i]]() ) ); @@ -485,7 +570,7 @@ void mesquiteMotionSolver::readOptions() ( new Mesquite::LPtoPTemplate ( - &qMetricTable_[qMetric_](), + &qMetricTable_[qMetric_[i]](), pValue, err ) @@ -501,7 +586,7 @@ void mesquiteMotionSolver::readOptions() ( new Mesquite::MaxTemplate ( - &qMetricTable_[qMetric_]() + &qMetricTable_[qMetric_[i]]() ) ); @@ -520,7 +605,7 @@ void mesquiteMotionSolver::readOptions() new Mesquite::PMeanPTemplate ( power, - &qMetricTable_[qMetric_]() + &qMetricTable_[qMetric_[i]]() ) ); @@ -533,7 +618,7 @@ void mesquiteMotionSolver::readOptions() ( new Mesquite::StdDevTemplate ( - &qMetricTable_[qMetric_]() + &qMetricTable_[qMetric_[i]]() ) ); @@ -546,7 +631,7 @@ void mesquiteMotionSolver::readOptions() ( new Mesquite::VarianceTemplate ( - &qMetricTable_[qMetric_]() + &qMetricTable_[qMetric_[i]]() ) ); @@ -564,7 +649,7 @@ void mesquiteMotionSolver::readOptions() new Mesquite::PatchPowerMeanP ( power, - &qMetricTable_[qMetric_]() + &qMetricTable_[qMetric_[i]]() ) ); @@ -711,7 +796,7 @@ void mesquiteMotionSolver::readOptions() ( dynamic_cast ( - &qMetricTable_[qMetric_]() + &qMetricTable_[qMetric_[0]]() ) ) ); @@ -1053,13 +1138,22 @@ void mesquiteMotionSolver::initArrays() return; } - // Prepare arrays for mesquite - vtxCoords_ = new double[3 * nPoints_]; - cellToNode_ = new unsigned long[4 * nCells_]; - fixFlags_ = new int[nPoints_]; + // Construct shape recognizers + hexMatcher hex; + tetMatcher tet; + prismMatcher prism; + pyrMatcher pyr; +// wedgeMatcher wedge; +// tetWedgeMatcher tetWedge; - // Set connectivity information - label cIndex = 0; + // Counters for different cell types + label nHex = 0; + label nTet = 0; + label nPrism = 0; + label nPyr = 0; +// label nWedge = 0; +// label nTetWedge = 0; + label nUnknown = 0; const faceList& meshFaces = mesh().faces(); const cellList& meshCells = mesh().cells(); @@ -1067,38 +1161,135 @@ void mesquiteMotionSolver::initArrays() forAll(meshCells, cellI) { - const cell& curCell = meshCells[cellI]; - const face& currFace = meshFaces[curCell[0]]; - const face& nextFace = meshFaces[curCell[1]]; - - // Get the fourth point - forAll(nextFace, pointI) + if (hex.isA(mesh(), cellI)) { - if - ( - nextFace[pointI] != currFace[0] - && nextFace[pointI] != currFace[1] - && nextFace[pointI] != currFace[2] - ) - { - // Fill in cellPoints in order - if (owner[curCell[0]] == cellI) - { - cellToNode_[cIndex++] = currFace[2]; - cellToNode_[cIndex++] = currFace[1]; - cellToNode_[cIndex++] = currFace[0]; - cellToNode_[cIndex++] = nextFace[pointI]; - } - else - { - cellToNode_[cIndex++] = currFace[0]; - cellToNode_[cIndex++] = currFace[1]; - cellToNode_[cIndex++] = currFace[2]; - cellToNode_[cIndex++] = nextFace[pointI]; - } + nHex++; + } + else if (tet.isA(mesh(), cellI)) + { + nTet++; + } + else if (pyr.isA(mesh(), cellI)) + { + nPyr++; + } + else if (prism.isA(mesh(), cellI)) + { + nPrism++; + } + // I think these are not supported by mesquite +// else if (wedge.isA(mesh(), cellI)) +// { +// nWedge++; +// } +// else if (tetWedge.isA(mesh(), cellI)) +// { +// nTetWedge++; +// } + else + { + // TODO: You will need to count all face points from these + // polyhedra to get the correct number of tets + // OR: improve the polyhedra support of mesquite +// nUnknown++; + Info << "You will need to count all face points from these " + << "polyhedra to get the correct number of tets!" + << endl; + } + } - break; + Info << "nHex = " << nHex << endl + << "nTet = " << nTet << endl + << "nPyr = " << nPyr << endl + << "nPrism = " << nPrism << endl + << "nUnknown = " << nUnknown << endl; + + // Prepare arrays for mesquite + vtxCoords_ = new double[3 * nPoints_]; + cellToNode_ = new unsigned long[8*nHex+4*nTet+5*nPyr+6*nPrism+nUnknown]; + fixFlags_ = new int[nPoints_]; + mixedTypes_ = new Mesquite::EntityTopology[nCells_]; + + // Set connectivity information + label cIndex = 0; + label cellIndex = 0; + + forAll(meshCells, cellI) + { + if (hex.isA(mesh(), cellI)) + { + hex.matchShape + ( + false, + meshFaces, + owner, + cellI, + meshCells[cellI] + ); + const labelList& pointLabelList = hex.vertLabels(); + forAll(pointLabelList, pointI) + { + cellToNode_[cIndex++] = pointLabelList[pointI]; } + mixedTypes_[cellIndex++] = Mesquite::HEXAHEDRON; + } + else if (tet.isA(mesh(), cellI)) + { + tet.matchShape + ( + false, + meshFaces, + owner, + cellI, + meshCells[cellI] + ); + const labelList& pointLabelList = tet.vertLabels(); + forAll(pointLabelList, pointI) + { + cellToNode_[cIndex++] = pointLabelList[pointI]; + } + mixedTypes_[cellIndex++] = Mesquite::TETRAHEDRON; + } + else if (pyr.isA(mesh(), cellI)) + { + pyr.matchShape + ( + false, + meshFaces, + owner, + cellI, + meshCells[cellI] + ); + // TODO: Check numbering of pyramids in mesquite! + const labelList& pointLabelList = pyr.vertLabels(); + forAll(pointLabelList, pointI) + { + cellToNode_[cIndex++] = pointLabelList[pointI]; + } + mixedTypes_[cellIndex++] = Mesquite::PYRAMID; + } + else if (prism.isA(mesh(), cellI)) + { + prism.matchShape + ( + false, + meshFaces, + owner, + cellI, + meshCells[cellI] + ); + const labelList& pointLabelList = prism.vertLabels(); + forAll(pointLabelList, pointI) + { + cellToNode_[cIndex++] = pointLabelList[pointI]; + } + mixedTypes_[cellIndex++] = Mesquite::PRISM; + } + else + { + // TODO: These cells have to be decomposed into tet4 elements + Info << "tet decomposition of polyhedra is not yet supported!" + << endl; } } @@ -1124,6 +1315,25 @@ void mesquiteMotionSolver::initArrays() } } + // Fetch the sub-dictionary + const dictionary& optionsDict = subDict("mesquiteOptions"); + + // fix additional points + if (optionsDict.found("fixPointsZone")) + { + word fixPoints = word(optionsDict.lookup("fixPointsZone")); + + label zoneID = mesh().pointZones().findZoneID(fixPoints); + if (zoneID > -1) + { + const pointZone& fixedPointLabels = mesh().pointZones()[zoneID]; + forAll(fixedPointLabels, pointI) + { + fixFlags_[fixedPointLabels[pointI]] = 1; + } + } + } + if (Pstream::parRun()) { initParallelConnectivity(); @@ -1901,6 +2111,21 @@ void mesquiteMotionSolver::solve() Mesquite::MsqError err; +// //- ArrayMesh object defined by Mesquite +// Mesquite::ArrayMesh +// msqMesh +// ( +// 3, // Number of coords per vertex +// nPoints_, // Number of vertices +// vtxCoords_, // The vertex coordinates +// fixFlags_, // Fixed vertex flags +// nCells_, // Number of elements +// Mesquite::TETRAHEDRON, // Element type +// cellToNode_, // Connectivity +// false, // Fortran-style array indexing +// 4 // Number of nodes per element +// ); + //- ArrayMesh object defined by Mesquite Mesquite::ArrayMesh msqMesh @@ -1910,10 +2135,10 @@ void mesquiteMotionSolver::solve() vtxCoords_, // The vertex coordinates fixFlags_, // Fixed vertex flags nCells_, // Number of elements - Mesquite::TETRAHEDRON, // Element type - cellToNode_, // Connectivity - false, // Fortran-style array indexing - 4 // Number of nodes per element + mixedTypes_, // array with all Element types + cellToNode_, // Connectivity - globel vertex indices + NULL, // Element offset connectivity + false // Fortran-style array indexing ); // Create an instruction queue @@ -1924,11 +2149,20 @@ void mesquiteMotionSolver::solve() optAlgorithm_->set_inner_termination_criterion(&tcInner_); // Set up the quality assessor - Mesquite::QualityAssessor qA(&qMetricTable_[qMetric_]()); + Mesquite::QualityAssessor qA(&qMetricTable_[qMetric_[0]]()); // Assess the quality of the initial mesh before smoothing queue.add_quality_assessor(&qA, err); + // Set up second optional quality assessor + Mesquite::QualityAssessor qASecond; + + if (qMetric_.size() == 2 ) + { + qASecond.add_quality_assessment(&qMetricTable_[qMetric_[1]]()); + queue.add_quality_assessor(&qASecond, err); + } + // Set the master quality improver queue.set_master_quality_improver ( @@ -1938,11 +2172,16 @@ void mesquiteMotionSolver::solve() // Assess the quality of the final mesh after smoothing queue.add_quality_assessor(&qA, err); + if (qMetric_.size() == 2 ) + { + queue.add_quality_assessor(&qASecond, err); + } // Disable slave output for parallel runs. if (Pstream::parRun() && !Pstream::master()) { qA.disable_printing_results(); + qASecond.disable_printing_results(); } // Launches optimization on the mesh diff --git a/src/dynamicMesh/meshMotion/mesquiteMotionSolver/mesquiteMotionSolver.H b/src/dynamicMesh/meshMotion/mesquiteMotionSolver/mesquiteMotionSolver.H index daeb63ac7..f8745cc5f 100644 --- a/src/dynamicMesh/meshMotion/mesquiteMotionSolver/mesquiteMotionSolver.H +++ b/src/dynamicMesh/meshMotion/mesquiteMotionSolver/mesquiteMotionSolver.H @@ -120,6 +120,9 @@ class mesquiteMotionSolver //- Flag array for vertices (fixed/free) mutable int* fixFlags_; + //- array with element types + mutable Mesquite::EntityTopology* mixedTypes_; + //- Reference points mutable pointIOField refPoints_; @@ -131,7 +134,7 @@ class mesquiteMotionSolver Map > auxSurfPointMap_; //- The quality metric - word qMetric_; + List qMetric_; //- Pointers to base element metric class HashTable > qMetricTable_; From 420ea04821a0c3d78510a5a9531eaa62a3a37749 Mon Sep 17 00:00:00 2001 From: "Oliver Borm (boroli)" Date: Tue, 3 May 2011 20:59:15 +0200 Subject: [PATCH 2/8] new version of mesquite solver, done by Sandeep Menon --- .../mesquiteMotionSolver.C | 2311 ++++++++++++++--- .../mesquiteMotionSolver.H | 91 +- 2 files changed, 2081 insertions(+), 321 deletions(-) diff --git a/src/dynamicMesh/meshMotion/mesquiteMotionSolver/mesquiteMotionSolver.C b/src/dynamicMesh/meshMotion/mesquiteMotionSolver/mesquiteMotionSolver.C index 40b3c0218..026bc5f5d 100644 --- a/src/dynamicMesh/meshMotion/mesquiteMotionSolver/mesquiteMotionSolver.C +++ b/src/dynamicMesh/meshMotion/mesquiteMotionSolver/mesquiteMotionSolver.C @@ -27,21 +27,22 @@ License #include "mesquiteMotionSolver.H" #include "Random.H" #include "IOmanip.H" +#include "SortableList.H" #include "globalMeshData.H" +#include "cyclicPolyPatch.H" #include "processorPolyPatch.H" #include "addToRunTimeSelectionTable.H" #include "polyMesh.H" +#include "matchPoints.H" #include "DimensionedField.H" #include "pointPatchField.H" #include "pointMesh.H" #include "hexMatcher.H" #include "tetMatcher.H" -#include "prismMatcher.H" #include "pyrMatcher.H" -// #include "wedgeMatcher.H" -// #include "tetWedgeMatcher.H" +#include "prismMatcher.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -63,12 +64,12 @@ mesquiteMotionSolver::mesquiteMotionSolver ) : motionSolver(mesh), + MeshObject(mesh), Mesh_(mesh), twoDMesh_(mesh.nGeometricD() == 2 ? true : false), + arraysInitialized_(false), nPoints_(mesh.nPoints()), nCells_(mesh.nCells()), - nAuxPoints_(0), - nAuxCells_(0), surfaceSmoothing_(false), volumeCorrection_(false), volCorrTolerance_(1e-20), @@ -77,10 +78,6 @@ mesquiteMotionSolver::mesquiteMotionSolver nSweeps_(1), surfInterval_(1), relax_(1.0), - vtxCoords_(NULL), - cellToNode_(NULL), - fixFlags_(NULL), - mixedTypes_(NULL), refPoints_ ( IOobject @@ -97,11 +94,9 @@ mesquiteMotionSolver::mesquiteMotionSolver { // Read options from the dictionary readOptions(); - - // Initialize connectivity arrays for Mesquite - initArrays(); } + mesquiteMotionSolver::mesquiteMotionSolver ( const polyMesh& mesh, @@ -109,12 +104,12 @@ mesquiteMotionSolver::mesquiteMotionSolver ) : motionSolver(mesh), + MeshObject(mesh), Mesh_(mesh), twoDMesh_(mesh.nGeometricD() == 2 ? true : false), + arraysInitialized_(false), nPoints_(mesh.nPoints()), nCells_(mesh.nCells()), - nAuxPoints_(0), - nAuxCells_(0), surfaceSmoothing_(false), volumeCorrection_(false), volCorrTolerance_(1e-20), @@ -123,10 +118,6 @@ mesquiteMotionSolver::mesquiteMotionSolver nSweeps_(1), surfInterval_(1), relax_(1.0), - vtxCoords_(NULL), - cellToNode_(NULL), - fixFlags_(NULL), - mixedTypes_(NULL), refPoints_ ( IOobject @@ -143,39 +134,94 @@ mesquiteMotionSolver::mesquiteMotionSolver { // Read options from the dictionary readOptions(); - - // Initialize connectivity arrays for Mesquite - initArrays(); } // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // mesquiteMotionSolver::~mesquiteMotionSolver() -{ - clearOut(); -} +{} -// Clear out addressing -void mesquiteMotionSolver::clearOut() -{ - // Delete memory pointers - delete [] vtxCoords_; - delete [] cellToNode_; - delete [] fixFlags_; - delete [] mixedTypes_; - - // Reset to NULL - vtxCoords_ = NULL; - cellToNode_ = NULL; - fixFlags_ = NULL; - mixedTypes_ = NULL; -} // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +// Parallel blocking send +void mesquiteMotionSolver::parWrite +( + const label toID, + const label& data +) +{ + OPstream::write + ( + Pstream::blocking, + toID, + reinterpret_cast(&data), + sizeof(label) + ); +} + + +// Parallel blocking receive +void mesquiteMotionSolver::parRead +( + const label fromID, + label& data +) +{ + IPstream::read + ( + Pstream::blocking, + fromID, + reinterpret_cast(&data), + sizeof(label) + ); +} + + +// Parallel non-blocking send for lists +template +void mesquiteMotionSolver::parWrite +( + const label toID, + const UList& data +) +{ + OPstream::write + ( + Pstream::nonBlocking, + toID, + reinterpret_cast(&data[0]), + data.size()*sizeof(Type) + ); +} + + +// Parallel non-blocking receive for lists +template +void mesquiteMotionSolver::parRead +( + const label fromID, + UList& data +) +{ + IPstream::read + ( + Pstream::nonBlocking, + fromID, + reinterpret_cast(&data[0]), + data.size()*sizeof(Type) + ); +} + + // Read options from the dictionary void mesquiteMotionSolver::readOptions() { + if (debug) + { + Info<< "Reading options for mesquiteMotionSolver" << endl; + } + // Fetch the sub-dictionary const dictionary& optionsDict = subDict("mesquiteOptions"); @@ -184,12 +230,21 @@ void mesquiteMotionSolver::readOptions() { labelHashSet slipPatchIDs; - // For 2D meshes, add all patches + const polyBoundaryMesh& boundary = mesh().boundaryMesh(); + + // For 2D meshes, add wedge / empty patches if (twoDMesh_) { - forAll(mesh().boundaryMesh(), patchI) + forAll(boundary, patchI) { - slipPatchIDs.insert(patchI); + if + ( + boundary[patchI].type() == "wedge" || + boundary[patchI].type() == "empty" + ) + { + slipPatchIDs.insert(patchI); + } } surfaceSmoothing_ = true; @@ -201,11 +256,7 @@ void mesquiteMotionSolver::readOptions() forAll(slipPatches, wordI) { word& patchName = slipPatches[wordI]; - - slipPatchIDs.insert - ( - mesh().boundaryMesh().findPatchID(patchName) - ); + slipPatchIDs.insert(boundary.findPatchID(patchName)); surfaceSmoothing_ = true; } @@ -267,8 +318,6 @@ void mesquiteMotionSolver::readOptions() optionsDict.subDict("coupledPatches") ); - const polyBoundaryMesh& boundary = mesh().boundaryMesh(); - // Determine master and slave patches forAllConstIter(dictionary, coupledPatches, dIter) { @@ -425,39 +474,64 @@ void mesquiteMotionSolver::readOptions() // ) // ); - qMetric_.setSize(1); + // Make sure we have at least one entry + qMetrics_.setSize(1); - if (optionsDict.found("optMetrics")) + // Read the optimization metric + if (optionsDict.isDict("optMetric")) { - const dictionary& optMetricsDict = optionsDict.subDict("optMetrics"); + // Metrics are in dictionary form + // does not work if two entries with the same name are present +// qMetrics_ = optionsDict.subDict("optMetric").toc(); - if (optMetricsDict.found("secondMetric")) + if (optionsDict.subDict("optMetric").found("secondMetric")) { - qMetric_.setSize(2); - qMetric_[1] = word(optMetricsDict.lookup("secondMetric")); + qMetrics_.setSize(2); + qMetrics_[1] = + word(optionsDict.subDict("optMetric").lookup("secondMetric")); } - if (optMetricsDict.found("firstMetric")) + qMetrics_[0] = + word(optionsDict.subDict("optMetric").lookup("firstMetric")); + + forAll(qMetrics_, wordI) { - qMetric_[0] = word(optMetricsDict.lookup("firstMetric")); + if (!qMetricTable_.found(qMetrics_[wordI])) + { + FatalErrorIn("void mesquiteMotionSolver::readOptions()") + << "Unrecognized quality metric: " + << qMetrics_[wordI] << nl + << "Available metrics are: " + << nl << qMetricTable_.toc() + << abort(FatalError); + } + else + { + Info<< "Selecting quality metric: " + << qMetrics_[wordI] << endl; + } } } else { - // Read the optimization metric - qMetric_[0] = word(optionsDict.lookup("optMetric")); - } + // Conventional keyword entry - if (!qMetricTable_.found(qMetric_[0])) - { - FatalErrorIn("void mesquiteMotionSolver::readOptions()") - << "Unrecognized quality metric: " << qMetric_[0] << nl - << "Available metrics are: " << nl << qMetricTable_.toc() - << abort(FatalError); - } - else - { - Info << "Selecting quality metric: " << qMetric_ << endl; + // Alias for convenience... + word& qMetric = qMetrics_[0]; + + qMetric = word(optionsDict.lookup("optMetric")); + + if (!qMetricTable_.found(qMetric)) + { + FatalErrorIn("void mesquiteMotionSolver::readOptions()") + << "Unrecognized quality metric: " << qMetric << nl + << "Available metrics are: " << nl << qMetricTable_.toc() + << abort(FatalError); + } + else + { + Info<< "Selecting quality metric: " << qMetric << endl; + } } // Define the objective function table, @@ -487,7 +561,7 @@ void mesquiteMotionSolver::readOptions() } else { - Info << "Selecting objective function: " << ofType << endl; + Info<< "Selecting objective function: " << ofType << endl; } // Instantiate the appropriate objective function @@ -512,6 +586,14 @@ void mesquiteMotionSolver::readOptions() << "Cannot make a composite of composite functions." << abort(FatalError); } + + // Make sure we have two quality metrics + if ( qMetrics_.size() < 2) + { + FatalErrorIn("void mesquiteMotionSolver::readOptions()") + << "Two quality metrics are needed for composite functions." + << abort(FatalError); + } } else // Check if a scaled function is requested @@ -554,7 +636,7 @@ void mesquiteMotionSolver::readOptions() ( new Mesquite::LInfTemplate ( - &qMetricTable_[qMetric_[i]]() + &qMetricTable_[qMetrics_[i]]() ) ); @@ -570,7 +652,7 @@ void mesquiteMotionSolver::readOptions() ( new Mesquite::LPtoPTemplate ( - &qMetricTable_[qMetric_[i]](), + &qMetricTable_[qMetrics_[i]](), pValue, err ) @@ -586,7 +668,7 @@ void mesquiteMotionSolver::readOptions() ( new Mesquite::MaxTemplate ( - &qMetricTable_[qMetric_[i]]() + &qMetricTable_[qMetrics_[i]]() ) ); @@ -605,7 +687,7 @@ void mesquiteMotionSolver::readOptions() new Mesquite::PMeanPTemplate ( power, - &qMetricTable_[qMetric_[i]]() + &qMetricTable_[qMetrics_[i]]() ) ); @@ -618,7 +700,7 @@ void mesquiteMotionSolver::readOptions() ( new Mesquite::StdDevTemplate ( - &qMetricTable_[qMetric_[i]]() + &qMetricTable_[qMetrics_[i]]() ) ); @@ -631,7 +713,7 @@ void mesquiteMotionSolver::readOptions() ( new Mesquite::VarianceTemplate ( - &qMetricTable_[qMetric_[i]]() + &qMetricTable_[qMetrics_[i]]() ) ); @@ -649,7 +731,7 @@ void mesquiteMotionSolver::readOptions() new Mesquite::PatchPowerMeanP ( power, - &qMetricTable_[qMetric_[i]]() + &qMetricTable_[qMetrics_[i]]() ) ); @@ -765,13 +847,13 @@ void mesquiteMotionSolver::readOptions() } else { - Info << "Selecting optimization algorithm: " << oaType << endl; + Info<< "Selecting optimization algorithm: " << oaType << endl; } // Instantiate the appropriate objective function label oaSelection = oaTable[oaType]; - switch (oaSelection) + switch(oaSelection) { case 0: { @@ -796,7 +878,7 @@ void mesquiteMotionSolver::readOptions() ( dynamic_cast ( - &qMetricTable_[qMetric_[0]]() + &qMetricTable_[qMetrics_[0]]() ) ) ); @@ -1014,9 +1096,9 @@ void mesquiteMotionSolver::readOptions() } else { - Info << "Inner termination criterion (tcInner) " - << "was not found. Using default values." - << endl; + Info<< "Inner termination criterion (tcInner) " + << "was not found. Using default values." + << endl; tcInner_.add_absolute_gradient_inf_norm(1e-4); } @@ -1063,9 +1145,9 @@ void mesquiteMotionSolver::readOptions() } else { - Info << "Outer termination criterion (tcOuter) " - << "was not found. Using default values." - << endl; + Info<< "Outer termination criterion (tcOuter) " + << "was not found. Using default values." + << endl; tcOuter_.add_iteration_limit(1); } @@ -1135,161 +1217,114 @@ void mesquiteMotionSolver::initArrays() if (twoDMesh_) { + // Initialise parallel connectivity, if necessary + initParallelConnectivity(); + + // Set the flag + arraysInitialized_ = true; + return; } // Construct shape recognizers - hexMatcher hex; - tetMatcher tet; - prismMatcher prism; - pyrMatcher pyr; -// wedgeMatcher wedge; -// tetWedgeMatcher tetWedge; - - // Counters for different cell types - label nHex = 0; - label nTet = 0; - label nPrism = 0; - label nPyr = 0; -// label nWedge = 0; -// label nTetWedge = 0; label nUnknown = 0; + FixedList nTypes(0); + FixedList, 4> matcher; + FixedList cellType; + + // Set types + matcher[0].set(new hexMatcher()); + cellType[0] = Mesquite::HEXAHEDRON; + + matcher[1].set(new tetMatcher()); + cellType[1] = Mesquite::TETRAHEDRON; + + matcher[2].set(new pyrMatcher()); + cellType[2] = Mesquite::PYRAMID; + + matcher[3].set(new prismMatcher()); + cellType[3] = Mesquite::PRISM; const faceList& meshFaces = mesh().faces(); const cellList& meshCells = mesh().cells(); - const labelList& owner = mesh().faceOwner(); + const labelList& faceOwner = mesh().faceOwner(); forAll(meshCells, cellI) { - if (hex.isA(mesh(), cellI)) + bool foundMatch = false; + + forAll(matcher, indexI) { - nHex++; + if (matcher[indexI]->isA(mesh(), cellI)) + { + nTypes[indexI]++; + + foundMatch = true; + break; + } } - else if (tet.isA(mesh(), cellI)) + + if (!foundMatch) { - nTet++; - } - else if (pyr.isA(mesh(), cellI)) - { - nPyr++; - } - else if (prism.isA(mesh(), cellI)) - { - nPrism++; - } - // I think these are not supported by mesquite -// else if (wedge.isA(mesh(), cellI)) -// { -// nWedge++; -// } -// else if (tetWedge.isA(mesh(), cellI)) -// { -// nTetWedge++; -// } - else - { - // TODO: You will need to count all face points from these - // polyhedra to get the correct number of tets - // OR: improve the polyhedra support of mesquite -// nUnknown++; - Info << "You will need to count all face points from these " - << "polyhedra to get the correct number of tets!" - << endl; + nUnknown++; } } - Info << "nHex = " << nHex << endl - << "nTet = " << nTet << endl - << "nPyr = " << nPyr << endl - << "nPrism = " << nPrism << endl - << "nUnknown = " << nUnknown << endl; + if (debug) + { + Pout<< "nHex = " << nTypes[0] << nl + << "nTet = " << nTypes[1] << nl + << "nPyr = " << nTypes[2] << nl + << "nPrism = " << nTypes[3] << nl + << "nUnknown = " << nUnknown << endl; + } + + label nCellToNode = + ( + (8 * nTypes[0]) + + (4 * nTypes[1]) + + (5 * nTypes[2]) + + (6 * nTypes[3]) + ); // Prepare arrays for mesquite - vtxCoords_ = new double[3 * nPoints_]; - cellToNode_ = new unsigned long[8*nHex+4*nTet+5*nPyr+6*nPrism+nUnknown]; - fixFlags_ = new int[nPoints_]; - mixedTypes_ = new Mesquite::EntityTopology[nCells_]; + vtxCoords_.setSize(3 * nPoints_); + cellToNode_.setSize(nCellToNode); + fixFlags_.setSize(nPoints_); + mixedTypes_.setSize(nCells_); // Set connectivity information - label cIndex = 0; - label cellIndex = 0; + label cIndex = 0, cellIndex = 0; forAll(meshCells, cellI) { - if (hex.isA(mesh(), cellI)) + forAll(matcher, indexI) { - hex.matchShape - ( - false, - meshFaces, - owner, - cellI, - meshCells[cellI] - ); - const labelList& pointLabelList = hex.vertLabels(); - forAll(pointLabelList, pointI) + if (matcher[indexI]->isA(mesh(), cellI)) { - cellToNode_[cIndex++] = pointLabelList[pointI]; + // Match shape + matcher[indexI]->matchShape + ( + false, + meshFaces, + faceOwner, + cellI, + meshCells[cellI] + ); + + // Set cellToNode + const labelList& cellToNode = matcher[indexI]->vertLabels(); + + forAll(cellToNode, nodeI) + { + cellToNode_[cIndex++] = cellToNode[nodeI]; + } + + // Set element type + mixedTypes_[cellIndex++] = cellType[indexI]; + + break; } - mixedTypes_[cellIndex++] = Mesquite::HEXAHEDRON; - } - else if (tet.isA(mesh(), cellI)) - { - tet.matchShape - ( - false, - meshFaces, - owner, - cellI, - meshCells[cellI] - ); - const labelList& pointLabelList = tet.vertLabels(); - forAll(pointLabelList, pointI) - { - cellToNode_[cIndex++] = pointLabelList[pointI]; - } - mixedTypes_[cellIndex++] = Mesquite::TETRAHEDRON; - } - else if (pyr.isA(mesh(), cellI)) - { - pyr.matchShape - ( - false, - meshFaces, - owner, - cellI, - meshCells[cellI] - ); - // TODO: Check numbering of pyramids in mesquite! - const labelList& pointLabelList = pyr.vertLabels(); - forAll(pointLabelList, pointI) - { - cellToNode_[cIndex++] = pointLabelList[pointI]; - } - mixedTypes_[cellIndex++] = Mesquite::PYRAMID; - } - else if (prism.isA(mesh(), cellI)) - { - prism.matchShape - ( - false, - meshFaces, - owner, - cellI, - meshCells[cellI] - ); - const labelList& pointLabelList = prism.vertLabels(); - forAll(pointLabelList, pointI) - { - cellToNode_[cIndex++] = pointLabelList[pointI]; - } - mixedTypes_[cellIndex++] = Mesquite::PRISM; - } - else - { - // TODO: These cells have to be decomposed into tet4 elements - Info << "tet decomposition of polyhedra is not yet supported!" - << endl; } } @@ -1315,43 +1350,1257 @@ void mesquiteMotionSolver::initArrays() } } - // Fetch the sub-dictionary + // Optionally fix additional zone points const dictionary& optionsDict = subDict("mesquiteOptions"); - // fix additional points if (optionsDict.found("fixPointsZone")) { - word fixPoints = word(optionsDict.lookup("fixPointsZone")); + wordList fixZones = optionsDict.subDict("fixPointsZone").toc(); - label zoneID = mesh().pointZones().findZoneID(fixPoints); - if (zoneID > -1) + forAll(fixZones, zoneI) { - const pointZone& fixedPointLabels = mesh().pointZones()[zoneID]; - forAll(fixedPointLabels, pointI) + label zoneID = mesh().pointZones().findZoneID(fixZones[zoneI]); + + if (zoneID > -1) { - fixFlags_[fixedPointLabels[pointI]] = 1; + const pointZone& pLabels = mesh().pointZones()[zoneID]; + + forAll(pLabels, pointI) + { + fixFlags_[pLabels[pointI]] = 1; + } } } } - if (Pstream::parRun()) + // Initialise parallel connectivity, if necessary + initParallelConnectivity(); + + // Set the flag + arraysInitialized_ = true; +} + + +// Identify coupled patches +void mesquiteMotionSolver::identifyCoupledPatches() +{ + if (!Pstream::parRun()) { - initParallelConnectivity(); + return; } + + if (debug) + { + Info<< "Identifying coupled patches" << endl; + } + + // Maintain a separate list of processor IDs in procIndices. + // This is done because this sub-domain may talk to processors + // that share only edges/points. + const polyBoundaryMesh& boundary = mesh().boundaryMesh(); + + // Fetch the list of global points from polyMesh. + const globalMeshData& gData = mesh().globalData(); + + const labelList& spA = gData.sharedPointAddr(); + const labelList& spL = gData.sharedPointLabels(); + + labelListList spB(Pstream::nProcs(), labelList(0)); + + if (gData.nGlobalPoints()) + { + if (debug) + { + Info<< " mesquiteMotionSolver::identifyCoupledPatches :" + << " Found " << gData.nGlobalPoints() + << " global points." << endl; + } + + // Send others my addressing. + for (label proc = 0; proc < Pstream::nProcs(); proc++) + { + if (proc != Pstream::myProcNo()) + { + // Send number of entities first. + parWrite(proc, spA.size()); + + // Send the buffer. + if (spA.size()) + { + parWrite(proc, spA); + } + } + } + + // Receive addressing from others + for (label proc = 0; proc < Pstream::nProcs(); proc++) + { + if (proc != Pstream::myProcNo()) + { + label procInfoSize = -1; + + // How many entities am I going to be receiving? + parRead(proc, procInfoSize); + + if (procInfoSize) + { + // Size the receive buffer. + spB[proc].setSize(procInfoSize, -1); + + // Schedule for receipt. + parRead(proc, spB[proc]); + } + } + } + } + else + if (debug) + { + Info<< " mesquiteMotionSolver::identifyCoupledPatches :" + << " Did not find any global points." << endl; + } + + labelHashSet immNeighbours; + labelListList procPoints(Pstream::nProcs()); + + // Track globally shared points + List > globalProcPoints + ( + Pstream::nProcs(), + DynamicList(5) + ); + + // Insert my immediate neighbours into the list. + forAll(boundary, pI) + { + if (isA(boundary[pI])) + { + const processorPolyPatch& pp = + ( + refCast(boundary[pI]) + ); + + label neiProcNo = pp.neighbProcNo(); + + // Insert all boundary points. + procPoints[neiProcNo] = pp.meshPoints(); + + // Keep track of immediate neighbours. + immNeighbours.insert(neiProcNo); + } + } + + if (gData.nGlobalPoints()) + { + // Wait for all transfers to complete. + OPstream::waitRequests(); + IPstream::waitRequests(); + + // Now loop through all processor addressing, and check if + // any labels coincide with my global shared points. + // If this is true, we need to be talking to that neighbour + // as well (if not already). + for (label proc = 0; proc < Pstream::nProcs(); proc++) + { + if (proc == Pstream::myProcNo()) + { + continue; + } + + bool foundGlobalMatch = false; + + // Fetch reference to buffer + const labelList& procBuffer = spB[proc]; + DynamicList& gPP = globalProcPoints[proc]; + + forAll(procBuffer, pointI) + { + forAll(spA, pointJ) + { + if (spA[pointJ] == procBuffer[pointI]) + { + // Make an entry, if one wasn't made already + if (findIndex(procPoints[proc], spL[pointJ]) == -1) + { + gPP.append(labelPair(spL[pointJ], spA[pointJ])); + } + + foundGlobalMatch = true; + + break; + } + } + } + + // Sort in ascending order of global point labels + if (gPP.size()) + { + SortableList