zhyang-dev/OpenFOAM-12
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

rigidBodyMeshMotion: Added support for mesh refinement/unrefinement with multiple bodies

Browse Source
This commit is contained in:
Henry Weller
2021-11-02 16:54:51 +00:00
parent a817efc9c6
commit a075e1a774
2 changed files with 147 additions and 114 deletions
Download Patch File Download Diff File Expand all files Collapse all files

251
src/rigidBodyMeshMotion/rigidBodyMeshMotion/rigidBodyMeshMotion.C
View File

@ -50,6 +50,57 @@ namespace Foam
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
Foam::List<Foam::septernion> Foam::rigidBodyMeshMotion::transforms0() const
{
List<septernion> 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::scalar>& Foam::rigidBodyMeshMotion::weights
(
const label pointi,
List<scalar>& 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<septernion> 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<septernion> transforms0(this->transforms0());
List<scalar> 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<const scalarField*> 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<septernion> transforms0(this->transforms0());
List<scalar> 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);

10
src/rigidBodyMeshMotion/rigidBodyMeshMotion/rigidBodyMeshMotion.H
View File

@ -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<septernion> 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<scalar>& weights(const label pointi, List<scalar>& w) const;
public: