diff --git a/src/rigidBodyMeshMotion/rigidBodyMeshMotion/rigidBodyMeshMotion.C b/src/rigidBodyMeshMotion/rigidBodyMeshMotion/rigidBodyMeshMotion.C index b1d2328630..ecabe04ff0 100644 --- a/src/rigidBodyMeshMotion/rigidBodyMeshMotion/rigidBodyMeshMotion.C +++ b/src/rigidBodyMeshMotion/rigidBodyMeshMotion/rigidBodyMeshMotion.C @@ -50,6 +50,57 @@ namespace Foam } +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +Foam::List Foam::rigidBodyMeshMotion::transforms0() const +{ + List transforms0(bodyMeshes_.size() + 1); + + forAll(bodyMeshes_, bi) + { + // Calculate the septernion equivalent of the transformation + transforms0[bi] = septernion(transform0(bodyMeshes_[bi].bodyID_)); + } + + transforms0[bodyMeshes_.size()] = septernion::I; + + return transforms0; +} + + +Foam::List& Foam::rigidBodyMeshMotion::weights +( + const label pointi, + List& w +) const +{ + // Initialise to 1 for the far-field weight + scalar sum1mw = 1; + + forAll(bodyMeshes_, bi) + { + w[bi] = bodyMeshes_[bi].weight_[pointi]; + sum1mw += w[bi]/(1 + small - w[bi]); + } + + // Calculate the limiter for wi/(1 - wi) to ensure the sum(wi) = 1 + scalar lambda = 1/sum1mw; + + // Limit wi/(1 - wi) and sum the resulting wi + scalar sumw = 0; + forAll(bodyMeshes_, bi) + { + w[bi] = lambda*w[bi]/(1 + small - w[bi]); + sumw += w[bi]; + } + + // Calculate the weight for the stationary far-field + w[bodyMeshes_.size()] = 1 - sumw; + + return w; +} + + // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // Foam::rigidBodyMeshMotion::bodyMesh::bodyMesh @@ -196,7 +247,10 @@ Type Foam::rigidBodyMeshMotion::bodyMesh::weight ) const { // Scaling: 1 up to di then linear down to 0 at do away from patches - Type weight = min(max((do_ - pDist)/(do_ - di_), scalar(0)), scalar(1)); + Type weight + ( + min(max((do_ - pDist)/(do_ - di_), scalar(0)), scalar(1)) + ); // Convert the weight function to a cosine weight = @@ -303,16 +357,13 @@ void Foam::rigidBodyMeshMotion::solve() } } + vectorField& pointDisplacement = pointDisplacement_.primitiveFieldRef(); + const pointField& points0 = this->points0(); + // Update the displacements if (bodyMeshes_.size() == 1) { - const septernion transform0 - ( - this->transform0(bodyMeshes_[0].bodyID_) - ); - - vectorField& pointDisplacement = pointDisplacement_.primitiveFieldRef(); - const pointField& points0 = this->points0(); + const septernion transform0(this->transform0(bodyMeshes_[0].bodyID_)); const scalarField& weight = bodyMeshes_[0].weight_; forAll(points0, pointi) @@ -340,48 +391,14 @@ void Foam::rigidBodyMeshMotion::solve() } else { - List transforms0(bodyMeshes_.size() + 1); - - forAll(bodyMeshes_, bi) - { - // Calculate the septernion equivalent of the transformation - transforms0[bi] = septernion(transform0(bodyMeshes_[bi].bodyID_)); - } - - transforms0[bodyMeshes_.size()] = septernion::I; - - vectorField& pointDisplacement = pointDisplacement_.primitiveFieldRef(); - const pointField& points0 = this->points0(); - + const List transforms0(this->transforms0()); List w(transforms0.size()); forAll(points0, pointi) { - // Initialise to 1 for the far-field weight - scalar sum1mw = 1; - - forAll(bodyMeshes_, bi) - { - w[bi] = bodyMeshes_[bi].weight_[pointi]; - sum1mw += w[bi]/(1 + small - w[bi]); - } - - // Calculate the limiter for wi/(1 - wi) to ensure the sum(wi) = 1 - scalar lambda = 1/sum1mw; - - // Limit wi/(1 - wi) and sum the resulting wi - scalar sumw = 0; - forAll(bodyMeshes_, bi) - { - w[bi] = lambda*w[bi]/(1 + small - w[bi]); - sumw += w[bi]; - } - - // Calculate the weight for the stationary far-field - w[bodyMeshes_.size()] = 1 - sumw; - pointDisplacement[pointi] = - average(transforms0, w).transformPoint(points0[pointi]) + average(transforms0, weights(pointi, w)) + .transformPoint(points0[pointi]) - points0[pointi]; } } @@ -409,16 +426,48 @@ void Foam::rigidBodyMeshMotion::updateMesh(const mapPolyMesh& mpm) ); const pointMesh& pMesh = pointMesh::New(mesh()); - pointField points0(points); + // Iterate to update the transformation of the new points to the + // corresponding points0, required because the body-point weights are + // calculated for points0 for (int iter=0; iter<3; iter++) { - // Calculate scaling factor everywhere for each meshed body - forAll(bodyMeshes_, bi) - { - const pointPatchDist pDist(pMesh, bodyMeshes_[bi].patchSet_, points0); - pointScalarField& weight = bodyMeshes_[bi].weight_; + // Calculate scaling factor everywhere for each meshed body + forAll(bodyMeshes_, bi) + { + const pointPatchDist pDist + ( + pMesh, + bodyMeshes_[bi].patchSet_, + points0 + ); + + pointScalarField& weight = bodyMeshes_[bi].weight_; + + forAll(points0, pointi) + { + const label oldPointi = mpm.pointMap()[pointi]; + + if (oldPointi >= 0) + { + if (mpm.reversePointMap()[oldPointi] != pointi) + { + weight[pointi] = bodyMeshes_[bi].weight(pDist[pointi]); + } + } + else + { + FatalErrorInFunction + << "Cannot determine co-ordinates " + "of introduced vertices." + << " New vertex " << pointi << " at co-ordinate " + << points[pointi] << exit(FatalError); + } + } + + pointConstraints::New(pMesh).constrain(weight); + } forAll(points0, pointi) { @@ -426,78 +475,52 @@ void Foam::rigidBodyMeshMotion::updateMesh(const mapPolyMesh& mpm) if (oldPointi >= 0) { - if (mpm.reversePointMap()[oldPointi] != pointi) + if (mpm.reversePointMap()[oldPointi] == pointi) { - weight[pointi] = bodyMeshes_[bi].weight(pDist[pointi]); - } - } - else - { - FatalErrorInFunction - << "Cannot determine co-ordinates of introduced vertices." - << " New vertex " << pointi << " at co-ordinate " - << points[pointi] << exit(FatalError); - } - } - - pointConstraints::New(pMesh).constrain(weight); - } - - forAll(points0, pointi) - { - const label oldPointi = mpm.pointMap()[pointi]; - - if (oldPointi >= 0) - { - if (mpm.reversePointMap()[oldPointi] == pointi) - { - points0[pointi] = points0_[oldPointi]; - } - else - { - if (bodyMeshes_.size() == 1) - { - // Use solid-body motion where weight = 1 - if (bodyMeshes_[0].weight_[pointi] > 1 - small) - { - points0[pointi] = - transform0(bodyMeshes_[0].bodyID_).inv() - .transformPoint(points[pointi]); - } - // Slerp septernion interpolation - else - { - points0[pointi] = - slerp - ( - septernion::I, - septernion(transform0(bodyMeshes_[0].bodyID_)), - bodyMeshes_[0].weight_[pointi] - ).invTransformPoint(points[pointi]); - } + points0[pointi] = points0_[oldPointi]; } else { - NotImplemented; - // labelList bodyIDs(bodyMeshes_.size()); - // List weights(bodyMeshes_.size()); - // forAll(bodyIDs, bi) - // { - // bodyIDs[bi] = bodyMeshes_[bi].bodyID_; - // weights[bi] = &bodyMeshes_[bi].weight_; - // } + if (bodyMeshes_.size() == 1) + { + // Use solid-body motion where weight = 1 + if (bodyMeshes_[0].weight_[pointi] > 1 - small) + { + points0[pointi] = + transform0(bodyMeshes_[0].bodyID_).inv() + .transformPoint(points[pointi]); + } + // Slerp septernion interpolation + else + { + points0[pointi] = + slerp + ( + septernion::I, + septernion + ( + transform0(bodyMeshes_[0].bodyID_) + ), + bodyMeshes_[0].weight_[pointi] + ).invTransformPoint(points[pointi]); + } + } + else + { + const List transforms0(this->transforms0()); + List w(transforms0.size()); - // points0[pointi] = invTransformPoints - // ( - // bodyIDs, - // weights, - // points[pointi] - // ); + forAll(points0, pointi) + { + points0[pointi] = + average(transforms0, weights(pointi, w)) + .invTransformPoint(points0[pointi]); + } + } } } } } - } points0_.transfer(points0); diff --git a/src/rigidBodyMeshMotion/rigidBodyMeshMotion/rigidBodyMeshMotion.H b/src/rigidBodyMeshMotion/rigidBodyMeshMotion/rigidBodyMeshMotion.H index 30ff295b0c..0f37514ee8 100644 --- a/src/rigidBodyMeshMotion/rigidBodyMeshMotion/rigidBodyMeshMotion.H +++ b/src/rigidBodyMeshMotion/rigidBodyMeshMotion/rigidBodyMeshMotion.H @@ -124,9 +124,19 @@ class rigidBodyMeshMotion //- Current time index (used for updating) label curTimeIndex_; + + // Private Member Functions + //- To avoid warning from clang using RBD::rigidBodyMotion::write; + //- Return the list of transformation for each body from time 0 + List transforms0() const; + + //- Return the list of weights for each body to pointi + // w is used as the work-space for the weights and returned + List& weights(const label pointi, List& w) const; + public: