COMP: resolve conflict
This commit is contained in:
andy
@ -184,6 +184,8 @@
|
||||
+ =setSet=: allows time range (e.g. 0:100) in combination with -batch argument
|
||||
to execute the commands for multiple times.
|
||||
* Post-processing
|
||||
+ =paraFoam=, =foamToVTK=: full support for polyhedral cell type in recent
|
||||
Paraview versions.
|
||||
+ =foamToEnsight=: parallel continuous data. new =-nodeValues= option to generate and output nodal
|
||||
field data.
|
||||
+ =singleCellMesh=: new utility to convert mesh and fields to a single cell
|
||||
|
||||
@ -218,7 +218,7 @@ void Foam::kineticTheoryModel::solve(const volTensorField& gradUat)
|
||||
// The solution is higly unstable close to the packing limit.
|
||||
gs0_ = radialModel_->g0
|
||||
(
|
||||
min(max(alpha_, 1e-6), alphaMax_ - 0.01),
|
||||
min(max(alpha_, scalar(1e-6)), alphaMax_ - 0.01),
|
||||
alphaMax_
|
||||
);
|
||||
|
||||
@ -255,7 +255,7 @@ void Foam::kineticTheoryModel::solve(const volTensorField& gradUat)
|
||||
volScalarField J1 = 3.0*betaPrim;
|
||||
volScalarField J2 =
|
||||
0.25*sqr(betaPrim)*da_*sqr(Ur)
|
||||
/(max(alpha_, 1e-6)*rhoa_*sqrtPi*(ThetaSqrt + TsmallSqrt));
|
||||
/(max(alpha_, scalar(1e-6))*rhoa_*sqrtPi*(ThetaSqrt + TsmallSqrt));
|
||||
|
||||
// bulk viscosity p. 45 (Lun et al. 1984).
|
||||
lambda_ = (4.0/3.0)*sqr(alpha_)*rhoa_*da_*gs0_*(1.0+e_)*ThetaSqrt/sqrtPi;
|
||||
@ -309,7 +309,11 @@ void Foam::kineticTheoryModel::solve(const volTensorField& gradUat)
|
||||
volScalarField t1 = K1*alpha_ + rhoa_;
|
||||
volScalarField l1 = -t1*trD;
|
||||
volScalarField l2 = sqr(t1)*tr2D;
|
||||
volScalarField l3 = 4.0*K4*max(alpha_, 1e-6)*(2.0*K3*trD2 + K2*tr2D);
|
||||
volScalarField l3 =
|
||||
4.0
|
||||
*K4
|
||||
*max(alpha_, scalar(1e-6))
|
||||
*(2.0*K3*trD2 + K2*tr2D);
|
||||
|
||||
Theta_ = sqr((l1 + sqrt(l2 + l3))/(2.0*(alpha_ + 1.0e-4)*K4));
|
||||
}
|
||||
|
||||
@ -1,5 +1,6 @@
|
||||
EXE_INC = \
|
||||
-I$(LIB_SRC)/meshTools/lnInclude
|
||||
-I$(LIB_SRC)/meshTools/lnInclude \
|
||||
-I$(LIB_SRC)/triSurface/lnInclude
|
||||
|
||||
EXE_LIBS = \
|
||||
-lmeshTools
|
||||
|
||||
@ -26,180 +26,38 @@ Application
|
||||
|
||||
Description
|
||||
Calculates the inertia tensor and principal axes and moments of a
|
||||
test face and tetrahedron.
|
||||
test face, tetrahedron and mesh.
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#include "argList.H"
|
||||
#include "Time.H"
|
||||
#include "polyMesh.H"
|
||||
#include "ListOps.H"
|
||||
#include "face.H"
|
||||
#include "tetPointRef.H"
|
||||
#include "triFaceList.H"
|
||||
#include "OFstream.H"
|
||||
#include "meshTools.H"
|
||||
#include "momentOfInertia.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
using namespace Foam;
|
||||
|
||||
void massPropertiesSolid
|
||||
(
|
||||
const pointField& pts,
|
||||
const triFaceList triFaces,
|
||||
scalar density,
|
||||
scalar& mass,
|
||||
vector& cM,
|
||||
tensor& J
|
||||
)
|
||||
{
|
||||
// Reimplemented from: Wm4PolyhedralMassProperties.cpp
|
||||
// File Version: 4.10.0 (2009/11/18)
|
||||
|
||||
// Geometric Tools, LC
|
||||
// Copyright (c) 1998-2010
|
||||
// Distributed under the Boost Software License, Version 1.0.
|
||||
// http://www.boost.org/LICENSE_1_0.txt
|
||||
// http://www.geometrictools.com/License/Boost/LICENSE_1_0.txt
|
||||
|
||||
// Boost Software License - Version 1.0 - August 17th, 2003
|
||||
|
||||
// Permission is hereby granted, free of charge, to any person or
|
||||
// organization obtaining a copy of the software and accompanying
|
||||
// documentation covered by this license (the "Software") to use,
|
||||
// reproduce, display, distribute, execute, and transmit the
|
||||
// Software, and to prepare derivative works of the Software, and
|
||||
// to permit third-parties to whom the Software is furnished to do
|
||||
// so, all subject to the following:
|
||||
|
||||
// The copyright notices in the Software and this entire
|
||||
// statement, including the above license grant, this restriction
|
||||
// and the following disclaimer, must be included in all copies of
|
||||
// the Software, in whole or in part, and all derivative works of
|
||||
// the Software, unless such copies or derivative works are solely
|
||||
// in the form of machine-executable object code generated by a
|
||||
// source language processor.
|
||||
|
||||
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
|
||||
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
|
||||
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE AND
|
||||
// NON-INFRINGEMENT. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR
|
||||
// ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE FOR ANY DAMAGES OR
|
||||
// OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
|
||||
// ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE
|
||||
// USE OR OTHER DEALINGS IN THE SOFTWARE.
|
||||
|
||||
const scalar r6 = 1.0/6.0;
|
||||
const scalar r24 = 1.0/24.0;
|
||||
const scalar r60 = 1.0/60.0;
|
||||
const scalar r120 = 1.0/120.0;
|
||||
|
||||
// order: 1, x, y, z, x^2, y^2, z^2, xy, yz, zx
|
||||
scalarField integrals(10, 0.0);
|
||||
|
||||
forAll(triFaces, i)
|
||||
{
|
||||
const triFace& tri(triFaces[i]);
|
||||
|
||||
// vertices of triangle i
|
||||
vector v0 = pts[tri[0]];
|
||||
vector v1 = pts[tri[1]];
|
||||
vector v2 = pts[tri[2]];
|
||||
|
||||
// cross product of edges
|
||||
vector eA = v1 - v0;
|
||||
vector eB = v2 - v0;
|
||||
vector n = eA ^ eB;
|
||||
|
||||
// compute integral terms
|
||||
scalar tmp0, tmp1, tmp2;
|
||||
|
||||
scalar f1x, f2x, f3x, g0x, g1x, g2x;
|
||||
|
||||
tmp0 = v0.x() + v1.x();
|
||||
f1x = tmp0 + v2.x();
|
||||
tmp1 = v0.x()*v0.x();
|
||||
tmp2 = tmp1 + v1.x()*tmp0;
|
||||
f2x = tmp2 + v2.x()*f1x;
|
||||
f3x = v0.x()*tmp1 + v1.x()*tmp2 + v2.x()*f2x;
|
||||
g0x = f2x + v0.x()*(f1x + v0.x());
|
||||
g1x = f2x + v1.x()*(f1x + v1.x());
|
||||
g2x = f2x + v2.x()*(f1x + v2.x());
|
||||
|
||||
scalar f1y, f2y, f3y, g0y, g1y, g2y;
|
||||
|
||||
tmp0 = v0.y() + v1.y();
|
||||
f1y = tmp0 + v2.y();
|
||||
tmp1 = v0.y()*v0.y();
|
||||
tmp2 = tmp1 + v1.y()*tmp0;
|
||||
f2y = tmp2 + v2.y()*f1y;
|
||||
f3y = v0.y()*tmp1 + v1.y()*tmp2 + v2.y()*f2y;
|
||||
g0y = f2y + v0.y()*(f1y + v0.y());
|
||||
g1y = f2y + v1.y()*(f1y + v1.y());
|
||||
g2y = f2y + v2.y()*(f1y + v2.y());
|
||||
|
||||
scalar f1z, f2z, f3z, g0z, g1z, g2z;
|
||||
|
||||
tmp0 = v0.z() + v1.z();
|
||||
f1z = tmp0 + v2.z();
|
||||
tmp1 = v0.z()*v0.z();
|
||||
tmp2 = tmp1 + v1.z()*tmp0;
|
||||
f2z = tmp2 + v2.z()*f1z;
|
||||
f3z = v0.z()*tmp1 + v1.z()*tmp2 + v2.z()*f2z;
|
||||
g0z = f2z + v0.z()*(f1z + v0.z());
|
||||
g1z = f2z + v1.z()*(f1z + v1.z());
|
||||
g2z = f2z + v2.z()*(f1z + v2.z());
|
||||
|
||||
// update integrals
|
||||
integrals[0] += n.x()*f1x;
|
||||
integrals[1] += n.x()*f2x;
|
||||
integrals[2] += n.y()*f2y;
|
||||
integrals[3] += n.z()*f2z;
|
||||
integrals[4] += n.x()*f3x;
|
||||
integrals[5] += n.y()*f3y;
|
||||
integrals[6] += n.z()*f3z;
|
||||
integrals[7] += n.x()*(v0.y()*g0x + v1.y()*g1x + v2.y()*g2x);
|
||||
integrals[8] += n.y()*(v0.z()*g0y + v1.z()*g1y + v2.z()*g2y);
|
||||
integrals[9] += n.z()*(v0.x()*g0z + v1.x()*g1z + v2.x()*g2z);
|
||||
}
|
||||
|
||||
integrals[0] *= r6;
|
||||
integrals[1] *= r24;
|
||||
integrals[2] *= r24;
|
||||
integrals[3] *= r24;
|
||||
integrals[4] *= r60;
|
||||
integrals[5] *= r60;
|
||||
integrals[6] *= r60;
|
||||
integrals[7] *= r120;
|
||||
integrals[8] *= r120;
|
||||
integrals[9] *= r120;
|
||||
|
||||
// mass
|
||||
mass = integrals[0];
|
||||
|
||||
// center of mass
|
||||
cM = vector(integrals[1], integrals[2], integrals[3])/mass;
|
||||
|
||||
// inertia relative to origin
|
||||
J.xx() = integrals[5] + integrals[6];
|
||||
J.xy() = -integrals[7];
|
||||
J.xz() = -integrals[9];
|
||||
J.yx() = J.xy();
|
||||
J.yy() = integrals[4] + integrals[6];
|
||||
J.yz() = -integrals[8];
|
||||
J.zx() = J.xz();
|
||||
J.zy() = J.yz();
|
||||
J.zz() = integrals[4] + integrals[5];
|
||||
|
||||
// inertia relative to center of mass
|
||||
J -= mass*((cM & cM)*I - cM*cM);
|
||||
|
||||
// Apply density
|
||||
mass *= density;
|
||||
J *= density;
|
||||
}
|
||||
|
||||
|
||||
int main(int argc, char *argv[])
|
||||
{
|
||||
argList::addOption
|
||||
(
|
||||
"cell",
|
||||
"label",
|
||||
"cell to use for inertia calculation, defaults to 0"
|
||||
);
|
||||
|
||||
#include "setRootCase.H"
|
||||
#include "createTime.H"
|
||||
#include "createPolyMesh.H"
|
||||
|
||||
scalar density = 1.0;
|
||||
|
||||
{
|
||||
@ -286,16 +144,7 @@ int main(int argc, char *argv[])
|
||||
vector cM = vector::zero;
|
||||
tensor J = tensor::zero;
|
||||
|
||||
massPropertiesSolid
|
||||
(
|
||||
|
||||
pts,
|
||||
tetFaces,
|
||||
density,
|
||||
m,
|
||||
cM,
|
||||
J
|
||||
);
|
||||
momentOfInertia::massPropertiesSolid(pts, tetFaces, density, m, cM, J);
|
||||
|
||||
vector eVal = eigenValues(J);
|
||||
|
||||
@ -344,7 +193,50 @@ int main(int argc, char *argv[])
|
||||
{
|
||||
str << "l " << nPts + 1 << ' ' << i + 1 << endl;
|
||||
}
|
||||
}
|
||||
|
||||
{
|
||||
const label cellI = args.optionLookupOrDefault("cell", 0);
|
||||
|
||||
tensorField mI = momentOfInertia::meshInertia(mesh);
|
||||
|
||||
tensor& J = mI[cellI];
|
||||
|
||||
vector eVal = eigenValues(J);
|
||||
|
||||
Info<< nl
|
||||
<< "Inertia tensor of cell " << cellI << " " << J << nl
|
||||
<< "eigenValues (principal moments) " << eVal << endl;
|
||||
|
||||
J /= cmptMax(eVal);
|
||||
|
||||
tensor eVec = eigenVectors(J);
|
||||
|
||||
Info<< "eigenVectors (principal axes, from normalised inertia) " << eVec
|
||||
<< endl;
|
||||
|
||||
OFstream str("cell_" + name(cellI) + "_inertia.obj");
|
||||
|
||||
Info<< nl << "Writing scaled principal axes of cell " << cellI << " to "
|
||||
<< str.name() << endl;
|
||||
|
||||
const point& cC = mesh.cellCentres()[cellI];
|
||||
|
||||
scalar scale = mag
|
||||
(
|
||||
(cC - mesh.faceCentres()[mesh.cells()[cellI][0]])
|
||||
/eVal.component(findMin(eVal))
|
||||
);
|
||||
|
||||
meshTools::writeOBJ(str, cC);
|
||||
meshTools::writeOBJ(str, cC + scale*eVal.x()*eVec.x());
|
||||
meshTools::writeOBJ(str, cC + scale*eVal.y()*eVec.y());
|
||||
meshTools::writeOBJ(str, cC + scale*eVal.z()*eVec.z());
|
||||
|
||||
for (label i = 1; i < 4; i++)
|
||||
{
|
||||
str << "l " << 1 << ' ' << i + 1 << endl;
|
||||
}
|
||||
}
|
||||
|
||||
Info<< nl << "End" << nl << endl;
|
||||
|
||||
@ -3,6 +3,7 @@ cd ${0%/*} || exit 1 # run from this directory
|
||||
set -x
|
||||
|
||||
wclean libso extrudeModel
|
||||
wclean
|
||||
wclean extrudeMesh
|
||||
wclean extrudeToRegionMesh
|
||||
|
||||
# ----------------------------------------------------------------- end-of-file
|
||||
@ -3,6 +3,8 @@ cd ${0%/*} || exit 1 # run from this directory
|
||||
set -x
|
||||
|
||||
wmake libso extrudeModel
|
||||
wmake
|
||||
wmake extrudeMesh
|
||||
wmake extrudeToRegionMesh
|
||||
|
||||
|
||||
# ----------------------------------------------------------------- end-of-file
|
||||
@ -1,6 +1,6 @@
|
||||
EXE_INC = \
|
||||
-IextrudedMesh \
|
||||
-IextrudeModel/lnInclude \
|
||||
-I../extrudeModel/lnInclude \
|
||||
-I$(LIB_SRC)/finiteVolume/lnInclude \
|
||||
-I$(LIB_SRC)/surfMesh/lnInclude \
|
||||
-I$(LIB_SRC)/meshTools/lnInclude \
|
||||
@ -214,9 +214,8 @@ int main(int argc, char *argv[])
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"extrudeProperties",
|
||||
runTimeExtruded.constant(),
|
||||
regionDir,
|
||||
"extrudeMeshDict",
|
||||
runTimeExtruded.system(),
|
||||
runTimeExtruded,
|
||||
IOobject::MUST_READ_IF_MODIFIED
|
||||
)
|
||||
@ -10,7 +10,7 @@ FoamFile
|
||||
version 2.0;
|
||||
format ascii;
|
||||
class dictionary;
|
||||
object extrudeProperties;
|
||||
object extrudeMeshDict;
|
||||
}
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
@ -1,9 +1,11 @@
|
||||
EXE_INC = \
|
||||
-I../extrudeModel/lnInclude \
|
||||
-I$(LIB_SRC)/finiteVolume/lnInclude \
|
||||
-I$(LIB_SRC)/meshTools/lnInclude \
|
||||
-I$(LIB_SRC)/dynamicMesh/lnInclude
|
||||
|
||||
EXE_LIBS = \
|
||||
-lextrudeModel \
|
||||
-lfiniteVolume \
|
||||
-lmeshTools \
|
||||
-ldynamicMesh
|
||||
@ -93,7 +93,7 @@ void Foam::createShellMesh::calcPointRegions
|
||||
label fp2 = findIndex(f2, pointI);
|
||||
label& region = pointRegions[face2][fp2];
|
||||
if (region != -1)
|
||||
{
|
||||
{
|
||||
FatalErrorIn
|
||||
(
|
||||
"createShellMesh::calcPointRegions(..)"
|
||||
@ -185,18 +185,20 @@ Foam::createShellMesh::createShellMesh
|
||||
|
||||
void Foam::createShellMesh::setRefinement
|
||||
(
|
||||
const pointField& thickness,
|
||||
const pointField& firstLayerDisp,
|
||||
const scalar expansionRatio,
|
||||
const label nLayers,
|
||||
const labelList& topPatchID,
|
||||
const labelList& bottomPatchID,
|
||||
const labelListList& extrudeEdgePatches,
|
||||
polyTopoChange& meshMod
|
||||
)
|
||||
{
|
||||
if (thickness.size() != regionPoints_.size())
|
||||
if (firstLayerDisp.size() != regionPoints_.size())
|
||||
{
|
||||
FatalErrorIn("createShellMesh::setRefinement(..)")
|
||||
<< "nRegions:" << regionPoints_.size()
|
||||
<< " thickness:" << thickness.size()
|
||||
<< " firstLayerDisp:" << firstLayerDisp.size()
|
||||
<< exit(FatalError);
|
||||
}
|
||||
|
||||
@ -224,30 +226,36 @@ void Foam::createShellMesh::setRefinement
|
||||
|
||||
|
||||
// From cell to patch (trivial)
|
||||
DynamicList<label> cellToFaceMap(patch_.size());
|
||||
DynamicList<label> cellToFaceMap(nLayers*patch_.size());
|
||||
// From face to patch+turning index
|
||||
DynamicList<label> faceToFaceMap(2*patch_.size()+patch_.nEdges());
|
||||
DynamicList<label> faceToFaceMap
|
||||
(
|
||||
(nLayers+1)*(patch_.size()+patch_.nEdges())
|
||||
);
|
||||
// From face to patch edge index
|
||||
DynamicList<label> faceToEdgeMap(patch_.nEdges()+patch_.nEdges());
|
||||
DynamicList<label> faceToEdgeMap(nLayers*(patch_.nEdges()+patch_.nEdges()));
|
||||
// From point to patch point index
|
||||
DynamicList<label> pointToPointMap(2*patch_.nPoints());
|
||||
DynamicList<label> pointToPointMap((nLayers+1)*patch_.nPoints());
|
||||
|
||||
|
||||
// Introduce new cell for every face
|
||||
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
labelList addedCells(patch_.size());
|
||||
labelList addedCells(nLayers*patch_.size());
|
||||
forAll(patch_, faceI)
|
||||
{
|
||||
addedCells[faceI] = meshMod.addCell
|
||||
(
|
||||
-1, // masterPointID
|
||||
-1, // masterEdgeID
|
||||
-1, // masterFaceID
|
||||
cellToFaceMap.size(), // masterCellID
|
||||
-1 // zoneID
|
||||
);
|
||||
cellToFaceMap.append(faceI);
|
||||
for (label layerI = 0; layerI < nLayers; layerI++)
|
||||
{
|
||||
addedCells[nLayers*faceI+layerI] = meshMod.addCell
|
||||
(
|
||||
-1, // masterPointID
|
||||
-1, // masterEdgeID
|
||||
-1, // masterFaceID
|
||||
cellToFaceMap.size(), // masterCellID
|
||||
-1 // zoneID
|
||||
);
|
||||
cellToFaceMap.append(faceI);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@ -261,7 +269,7 @@ void Foam::createShellMesh::setRefinement
|
||||
meshMod.addPoint
|
||||
(
|
||||
patch_.localPoints()[pointI], // point
|
||||
pointToPointMap.size(), // masterPointID
|
||||
pointToPointMap.size(), // masterPointID
|
||||
-1, // zoneID
|
||||
true // inCell
|
||||
);
|
||||
@ -277,25 +285,28 @@ void Foam::createShellMesh::setRefinement
|
||||
// Introduce new points (one for every region)
|
||||
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
labelList addedPoints(regionPoints_.size());
|
||||
labelList addedPoints(nLayers*regionPoints_.size());
|
||||
forAll(regionPoints_, regionI)
|
||||
{
|
||||
label pointI = regionPoints_[regionI];
|
||||
point extrudedPt = patch_.localPoints()[pointI] + thickness[regionI];
|
||||
|
||||
addedPoints[regionI] = meshMod.addPoint
|
||||
(
|
||||
extrudedPt, // point
|
||||
pointToPointMap.size(), // masterPointID - used only addressing
|
||||
-1, // zoneID
|
||||
true // inCell
|
||||
);
|
||||
pointToPointMap.append(pointI);
|
||||
point pt = patch_.localPoints()[pointI];
|
||||
point disp = firstLayerDisp[regionI];
|
||||
for (label layerI = 0; layerI < nLayers; layerI++)
|
||||
{
|
||||
pt += disp;
|
||||
|
||||
//Pout<< "Added top point " << addedPoints[regionI]
|
||||
// << " at " << extrudedPt
|
||||
// << " from point " << pointI
|
||||
// << endl;
|
||||
addedPoints[nLayers*regionI+layerI] = meshMod.addPoint
|
||||
(
|
||||
pt, // point
|
||||
pointToPointMap.size(), // masterPointID - used only addressing
|
||||
-1, // zoneID
|
||||
true // inCell
|
||||
);
|
||||
pointToPointMap.append(pointI);
|
||||
|
||||
disp *= expansionRatio;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@ -305,61 +316,84 @@ void Foam::createShellMesh::setRefinement
|
||||
meshMod.addFace
|
||||
(
|
||||
patch_.localFaces()[faceI].reverseFace(),// vertices
|
||||
addedCells[faceI], // own
|
||||
addedCells[nLayers*faceI], // own
|
||||
-1, // nei
|
||||
-1, // masterPointID
|
||||
-1, // masterEdgeID
|
||||
faceToFaceMap.size(), // masterFaceID : current faceI
|
||||
true, // flipFaceFlux
|
||||
bottomPatchID[faceI],// patchID
|
||||
bottomPatchID[faceI], // patchID
|
||||
-1, // zoneID
|
||||
false // zoneFlip
|
||||
);
|
||||
faceToFaceMap.append(-faceI-1); // points to flipped original face
|
||||
faceToEdgeMap.append(-1);
|
||||
|
||||
//const face newF(patch_.localFaces()[faceI].reverseFace());
|
||||
//Pout<< "Added bottom face "
|
||||
// << patch_.localFaces()[faceI].reverseFace()
|
||||
// << newF
|
||||
// << " coords:" << UIndirectList<point>(meshMod.points(), newF)
|
||||
// << " own " << addedCells[faceI]
|
||||
// << " patch:" << bottomPatchID[faceI]
|
||||
// << " at " << patch_.faceCentres()[faceI]
|
||||
// << endl;
|
||||
|
||||
}
|
||||
|
||||
// Add face on top
|
||||
// Add inbetween faces and face on top
|
||||
forAll(patch_.localFaces(), faceI)
|
||||
{
|
||||
// Get face in original ordering
|
||||
const face& f = patch_.localFaces()[faceI];
|
||||
|
||||
// Pick up point based on region
|
||||
face newF(f.size());
|
||||
forAll(f, fp)
|
||||
|
||||
for (label layerI = 0; layerI < nLayers; layerI++)
|
||||
{
|
||||
label region = pointRegions_[faceI][fp];
|
||||
newF[fp] = addedPoints[region];
|
||||
// Pick up point based on region and layer
|
||||
forAll(f, fp)
|
||||
{
|
||||
label region = pointRegions_[faceI][fp];
|
||||
newF[fp] = addedPoints[region*nLayers+layerI];
|
||||
}
|
||||
|
||||
label own = addedCells[faceI*nLayers+layerI];
|
||||
label nei;
|
||||
label patchI;
|
||||
if (layerI == nLayers-1)
|
||||
{
|
||||
nei = -1;
|
||||
patchI = topPatchID[faceI];
|
||||
}
|
||||
else
|
||||
{
|
||||
nei = addedCells[faceI*nLayers+layerI+1];
|
||||
patchI = -1;
|
||||
}
|
||||
|
||||
meshMod.addFace
|
||||
(
|
||||
newF, // vertices
|
||||
own, // own
|
||||
nei, // nei
|
||||
-1, // masterPointID
|
||||
-1, // masterEdgeID
|
||||
faceToFaceMap.size(), // masterFaceID : current faceI
|
||||
false, // flipFaceFlux
|
||||
patchI, // patchID
|
||||
-1, // zoneID
|
||||
false // zoneFlip
|
||||
);
|
||||
faceToFaceMap.append(faceI+1); // unflipped
|
||||
faceToEdgeMap.append(-1);
|
||||
|
||||
//Pout<< "Added inbetween face " << newF
|
||||
// << " coords:" << UIndirectList<point>(meshMod.points(), newF)
|
||||
// << " at layer " << layerI
|
||||
// << " own " << own
|
||||
// << " nei " << nei
|
||||
// << " at " << patch_.faceCentres()[faceI]
|
||||
// << endl;
|
||||
}
|
||||
|
||||
meshMod.addFace
|
||||
(
|
||||
newF, // vertices
|
||||
addedCells[faceI], // own
|
||||
-1, // nei
|
||||
-1, // masterPointID
|
||||
-1, // masterEdgeID
|
||||
faceToFaceMap.size(), // masterFaceID : current faceI
|
||||
false, // flipFaceFlux
|
||||
topPatchID[faceI], // patchID
|
||||
-1, // zoneID
|
||||
false // zoneFlip
|
||||
);
|
||||
faceToFaceMap.append(faceI+1); // unflipped
|
||||
faceToEdgeMap.append(-1);
|
||||
|
||||
//Pout<< "Added top face " << newF
|
||||
// << " own " << addedCells[faceI]
|
||||
// << " at " << patch_.faceCentres()[faceI]
|
||||
// << endl;
|
||||
}
|
||||
|
||||
|
||||
@ -376,8 +410,7 @@ void Foam::createShellMesh::setRefinement
|
||||
|
||||
if (ePatches.size() == 0)
|
||||
{
|
||||
// internal face.
|
||||
|
||||
// Internal face
|
||||
if (eFaces.size() != 2)
|
||||
{
|
||||
FatalErrorIn("createShellMesh::setRefinement(..)")
|
||||
@ -385,61 +418,6 @@ void Foam::createShellMesh::setRefinement
|
||||
<< " not internal but does not have side-patches defined."
|
||||
<< exit(FatalError);
|
||||
}
|
||||
|
||||
// Extrude
|
||||
|
||||
// Make face pointing in to eFaces[0] so out of new master face
|
||||
const face& f = patch_.localFaces()[eFaces[0]];
|
||||
const edge& e = patch_.edges()[edgeI];
|
||||
|
||||
label fp0 = findIndex(f, e[0]);
|
||||
label fp1 = f.fcIndex(fp0);
|
||||
|
||||
if (f[fp1] != e[1])
|
||||
{
|
||||
fp1 = fp0;
|
||||
fp0 = f.rcIndex(fp1);
|
||||
}
|
||||
|
||||
face newF(4);
|
||||
newF[0] = f[fp0];
|
||||
newF[1] = f[fp1];
|
||||
newF[2] = addedPoints[pointRegions_[eFaces[0]][fp1]];
|
||||
newF[3] = addedPoints[pointRegions_[eFaces[0]][fp0]];
|
||||
|
||||
label minCellI = addedCells[eFaces[0]];
|
||||
label maxCellI = addedCells[eFaces[1]];
|
||||
|
||||
if (minCellI > maxCellI)
|
||||
{
|
||||
// Swap
|
||||
Swap(minCellI, maxCellI);
|
||||
newF = newF.reverseFace();
|
||||
}
|
||||
|
||||
//Pout<< "for internal edge:" << e
|
||||
// << " at:" << patch_.localPoints()[e[0]]
|
||||
// << patch_.localPoints()[e[1]]
|
||||
// << " adding face:" << newF
|
||||
// << " from f:" << f
|
||||
// << " inbetween " << minCellI << " and " << maxCellI << endl;
|
||||
|
||||
// newF already outwards pointing.
|
||||
meshMod.addFace
|
||||
(
|
||||
newF, // vertices
|
||||
minCellI, // own
|
||||
maxCellI, // nei
|
||||
-1, // masterPointID
|
||||
-1, // masterEdgeID
|
||||
faceToFaceMap.size(), // masterFaceID
|
||||
false, // flipFaceFlux
|
||||
-1, // patchID
|
||||
-1, // zoneID
|
||||
false // zoneFlip
|
||||
);
|
||||
faceToFaceMap.append(0);
|
||||
faceToEdgeMap.append(edgeI);
|
||||
}
|
||||
else
|
||||
{
|
||||
@ -451,52 +429,91 @@ void Foam::createShellMesh::setRefinement
|
||||
<< " but only " << ePatches.size()
|
||||
<< " boundary faces defined." << exit(FatalError);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// Extrude eFaces[0]
|
||||
label minFaceI = eFaces[0];
|
||||
|
||||
// Make face pointing in to eFaces[0] so out of new master face
|
||||
const face& f = patch_.localFaces()[minFaceI];
|
||||
// Make face pointing in to eFaces[0] so out of new master face
|
||||
const face& f = patch_.localFaces()[eFaces[0]];
|
||||
const edge& e = patch_.edges()[edgeI];
|
||||
|
||||
const edge& e = patch_.edges()[edgeI];
|
||||
label fp0 = findIndex(f, e[0]);
|
||||
label fp1 = f.fcIndex(fp0);
|
||||
label fp0 = findIndex(f, e[0]);
|
||||
label fp1 = f.fcIndex(fp0);
|
||||
|
||||
if (f[fp1] != e[1])
|
||||
if (f[fp1] != e[1])
|
||||
{
|
||||
fp1 = fp0;
|
||||
fp0 = f.rcIndex(fp1);
|
||||
}
|
||||
|
||||
face newF(4);
|
||||
|
||||
for (label layerI = 0; layerI < nLayers; layerI++)
|
||||
{
|
||||
label region0 = pointRegions_[eFaces[0]][fp0];
|
||||
label region1 = pointRegions_[eFaces[0]][fp1];
|
||||
|
||||
if (layerI == 0)
|
||||
{
|
||||
fp1 = fp0;
|
||||
fp0 = f.rcIndex(fp1);
|
||||
newF[0] = f[fp0];
|
||||
newF[1] = f[fp1];
|
||||
newF[2] = addedPoints[nLayers*region1+layerI];
|
||||
newF[3] = addedPoints[nLayers*region0+layerI];
|
||||
}
|
||||
else
|
||||
{
|
||||
newF[0] = addedPoints[nLayers*region0+layerI-1];
|
||||
newF[1] = addedPoints[nLayers*region1+layerI-1];
|
||||
newF[2] = addedPoints[nLayers*region1+layerI];
|
||||
newF[3] = addedPoints[nLayers*region0+layerI];
|
||||
}
|
||||
|
||||
face newF(4);
|
||||
newF[0] = f[fp0];
|
||||
newF[1] = f[fp1];
|
||||
newF[2] = addedPoints[pointRegions_[minFaceI][fp1]];
|
||||
newF[3] = addedPoints[pointRegions_[minFaceI][fp0]];
|
||||
label minCellI = addedCells[nLayers*eFaces[0]+layerI];
|
||||
label maxCellI;
|
||||
label patchI;
|
||||
if (ePatches.size() == 0)
|
||||
{
|
||||
maxCellI = addedCells[nLayers*eFaces[1]+layerI];
|
||||
if (minCellI > maxCellI)
|
||||
{
|
||||
// Swap
|
||||
Swap(minCellI, maxCellI);
|
||||
newF = newF.reverseFace();
|
||||
}
|
||||
patchI = -1;
|
||||
}
|
||||
else
|
||||
{
|
||||
maxCellI = -1;
|
||||
patchI = ePatches[0];
|
||||
}
|
||||
|
||||
|
||||
//Pout<< "for external edge:" << e
|
||||
// << " at:" << patch_.localPoints()[e[0]]
|
||||
// << patch_.localPoints()[e[1]]
|
||||
// << " adding first patch face:" << newF
|
||||
// << " from:" << f
|
||||
// << " into patch:" << ePatches[0]
|
||||
// << " own:" << addedCells[minFaceI]
|
||||
// << endl;
|
||||
//{
|
||||
// Pout<< "Adding from face:" << patch_.faceCentres()[eFaces[0]]
|
||||
// << " from edge:"
|
||||
// << patch_.localPoints()[f[fp0]]
|
||||
// << patch_.localPoints()[f[fp1]]
|
||||
// << " at layer:" << layerI
|
||||
// << " with new points:" << newF
|
||||
// << " locations:"
|
||||
// << UIndirectList<point>(meshMod.points(), newF)
|
||||
// << " own:" << minCellI
|
||||
// << " nei:" << maxCellI
|
||||
// << endl;
|
||||
//}
|
||||
|
||||
|
||||
// newF already outwards pointing.
|
||||
meshMod.addFace
|
||||
(
|
||||
newF, // vertices
|
||||
addedCells[minFaceI], // own
|
||||
-1, // nei
|
||||
minCellI, // own
|
||||
maxCellI, // nei
|
||||
-1, // masterPointID
|
||||
-1, // masterEdgeID
|
||||
faceToFaceMap.size(), // masterFaceID
|
||||
false, // flipFaceFlux
|
||||
ePatches[0], // patchID
|
||||
patchI, // patchID
|
||||
-1, // zoneID
|
||||
false // zoneFlip
|
||||
);
|
||||
@ -532,35 +549,57 @@ void Foam::createShellMesh::setRefinement
|
||||
}
|
||||
|
||||
face newF(4);
|
||||
newF[0] = f[fp0];
|
||||
newF[1] = f[fp1];
|
||||
newF[2] = addedPoints[pointRegions_[minFaceI][fp1]];
|
||||
newF[3] = addedPoints[pointRegions_[minFaceI][fp0]];
|
||||
for (label layerI = 0; layerI < nLayers; layerI++)
|
||||
{
|
||||
label region0 = pointRegions_[minFaceI][fp0];
|
||||
label region1 = pointRegions_[minFaceI][fp1];
|
||||
|
||||
//Pout<< "for external edge:" << e
|
||||
// << " at:" << patch_.localPoints()[e[0]]
|
||||
// << patch_.localPoints()[e[1]]
|
||||
// << " adding patch face:" << newF
|
||||
// << " from:" << f
|
||||
// << " into patch:" << ePatches[i]
|
||||
// << endl;
|
||||
if (layerI == 0)
|
||||
{
|
||||
newF[0] = f[fp0];
|
||||
newF[1] = f[fp1];
|
||||
newF[2] = addedPoints[nLayers*region1+layerI];
|
||||
newF[3] = addedPoints[nLayers*region0+layerI];
|
||||
}
|
||||
else
|
||||
{
|
||||
newF[0] = addedPoints[nLayers*region0+layerI-1];
|
||||
newF[1] = addedPoints[nLayers*region1+layerI-1];
|
||||
newF[2] = addedPoints[nLayers*region1+layerI];
|
||||
newF[3] = addedPoints[nLayers*region0+layerI];
|
||||
}
|
||||
|
||||
// newF already outwards pointing.
|
||||
meshMod.addFace
|
||||
(
|
||||
newF, // vertices
|
||||
addedCells[minFaceI], // own
|
||||
-1, // nei
|
||||
-1, // masterPointID
|
||||
-1, // masterEdgeID
|
||||
faceToFaceMap.size(), // masterFaceID
|
||||
false, // flipFaceFlux
|
||||
ePatches[i], // patchID
|
||||
-1, // zoneID
|
||||
false // zoneFlip
|
||||
);
|
||||
faceToFaceMap.append(0);
|
||||
faceToEdgeMap.append(edgeI);
|
||||
////if (ePatches.size() == 0)
|
||||
//{
|
||||
// Pout<< "Adding from MULTI face:"
|
||||
// << patch_.faceCentres()[minFaceI]
|
||||
// << " from edge:"
|
||||
// << patch_.localPoints()[f[fp0]]
|
||||
// << patch_.localPoints()[f[fp1]]
|
||||
// << " at layer:" << layerI
|
||||
// << " with new points:" << newF
|
||||
// << " locations:"
|
||||
// << UIndirectList<point>(meshMod.points(), newF)
|
||||
// << endl;
|
||||
//}
|
||||
|
||||
// newF already outwards pointing.
|
||||
meshMod.addFace
|
||||
(
|
||||
newF, // vertices
|
||||
addedCells[nLayers*minFaceI+layerI], // own
|
||||
-1, // nei
|
||||
-1, // masterPointID
|
||||
-1, // masterEdgeID
|
||||
faceToFaceMap.size(), // masterFaceID
|
||||
false, // flipFaceFlux
|
||||
ePatches[i], // patchID
|
||||
-1, // zoneID
|
||||
false // zoneFlip
|
||||
);
|
||||
faceToFaceMap.append(0);
|
||||
faceToEdgeMap.append(edgeI);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -109,7 +109,8 @@ public:
|
||||
|
||||
// Access
|
||||
|
||||
//- From region cell to patch face
|
||||
//- From region cell to patch face. Consecutively added so
|
||||
// cell at layerI is at patchFaceI*nLayers+layerI
|
||||
const labelList& cellToFaceMap() const
|
||||
{
|
||||
return cellToFaceMap_;
|
||||
@ -120,7 +121,7 @@ public:
|
||||
// be in top patch
|
||||
// < 0 : face in opposite orientation as patch face. face will
|
||||
// be in bottom patch
|
||||
// = 0 : for all side faces
|
||||
// = 0 : for all side and internal faces
|
||||
const labelList& faceToFaceMap() const
|
||||
{
|
||||
return faceToFaceMap_;
|
||||
@ -153,7 +154,9 @@ public:
|
||||
//- Play commands into polyTopoChange to create layer mesh.
|
||||
void setRefinement
|
||||
(
|
||||
const pointField& thickness,
|
||||
const pointField& firstLayerThickness,
|
||||
const scalar expansionRatio,
|
||||
const label nLayers,
|
||||
const labelList& topPatchID,
|
||||
const labelList& bottomPatchID,
|
||||
const labelListList& extrudeEdgePatches,
|
||||
@ -133,6 +133,7 @@ Usage
|
||||
#include "syncTools.H"
|
||||
#include "cyclicPolyPatch.H"
|
||||
#include "nonuniformTransformCyclicPolyPatch.H"
|
||||
#include "extrudeModel.H"
|
||||
|
||||
using namespace Foam;
|
||||
|
||||
@ -689,39 +690,6 @@ void countExtrudePatches
|
||||
}
|
||||
|
||||
|
||||
// Lexical ordering for vectors.
|
||||
bool lessThan(const point& x, const point& y)
|
||||
{
|
||||
for (direction dir = 0; dir < point::nComponents; dir++)
|
||||
{
|
||||
if (x[dir] < y[dir]) return true;
|
||||
if (x[dir] > y[dir]) return false;
|
||||
}
|
||||
return false;
|
||||
}
|
||||
|
||||
|
||||
// Combine vectors
|
||||
class minEqVectorOp
|
||||
{
|
||||
public:
|
||||
void operator()(vector& x, const vector& y) const
|
||||
{
|
||||
if (y != vector::zero)
|
||||
{
|
||||
if (x == vector::zero)
|
||||
{
|
||||
x = y;
|
||||
}
|
||||
else if (lessThan(y, x))
|
||||
{
|
||||
x = y;
|
||||
}
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
// Constrain&sync normals on points that are on coupled patches to make sure
|
||||
// the face extruded from the edge has a valid normal with its coupled
|
||||
// equivalent.
|
||||
@ -846,14 +814,14 @@ void constrainCoupledNormals
|
||||
(
|
||||
mesh,
|
||||
edgeNormals0,
|
||||
minEqVectorOp(),
|
||||
maxMagSqrEqOp<vector>(),
|
||||
vector::zero // nullValue
|
||||
);
|
||||
syncTools::syncEdgeList
|
||||
(
|
||||
mesh,
|
||||
edgeNormals1,
|
||||
minEqVectorOp(),
|
||||
maxMagSqrEqOp<vector>(),
|
||||
vector::zero // nullValue
|
||||
);
|
||||
|
||||
@ -951,14 +919,9 @@ tmp<pointField> calcOffset
|
||||
int main(int argc, char *argv[])
|
||||
{
|
||||
argList::noParallel();
|
||||
argList::validArgs.append("shellRegion");
|
||||
argList::validArgs.append("faceZones");
|
||||
argList::validArgs.append("thickness");
|
||||
|
||||
Foam::argList::addBoolOption
|
||||
argList::addNote
|
||||
(
|
||||
"oneD",
|
||||
"generate columns of 1D cells"
|
||||
"Create region mesh by extruding a faceZone"
|
||||
);
|
||||
|
||||
#include "addRegionOption.H"
|
||||
@ -966,19 +929,34 @@ int main(int argc, char *argv[])
|
||||
#include "setRootCase.H"
|
||||
#include "createTime.H"
|
||||
#include "createNamedMesh.H"
|
||||
const word oldInstance = mesh.pointsInstance();
|
||||
|
||||
word shellRegionName = args.additionalArgs()[0];
|
||||
const wordList zoneNames(IStringStream(args.additionalArgs()[1])());
|
||||
scalar thickness = readScalar(IStringStream(args.additionalArgs()[2])());
|
||||
const word oldInstance = mesh.pointsInstance();
|
||||
bool overwrite = args.optionFound("overwrite");
|
||||
bool oneD = args.optionFound("oneD");
|
||||
|
||||
IOdictionary dict
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"extrudeToRegionMeshDict",
|
||||
runTime.system(),
|
||||
runTime,
|
||||
IOobject::MUST_READ_IF_MODIFIED
|
||||
)
|
||||
);
|
||||
|
||||
// Point generator
|
||||
autoPtr<extrudeModel> model(extrudeModel::New(dict));
|
||||
|
||||
|
||||
// Region
|
||||
const word shellRegionName(dict.lookup("region"));
|
||||
const wordList zoneNames(dict.lookup("faceZones"));
|
||||
const Switch oneD(dict.lookup("oneD"));
|
||||
|
||||
|
||||
Info<< "Extruding zones " << zoneNames
|
||||
<< " on mesh " << regionName
|
||||
<< " into shell mesh " << shellRegionName
|
||||
<< " of thickness " << thickness << nl
|
||||
<< endl;
|
||||
|
||||
if (shellRegionName == regionName)
|
||||
@ -1081,7 +1059,6 @@ int main(int argc, char *argv[])
|
||||
);
|
||||
|
||||
|
||||
|
||||
// Check whether the zone is internal or external faces to determine
|
||||
// what patch type to insert. Cannot be mixed
|
||||
// since then how to couple? - directMapped only valid for a whole patch.
|
||||
@ -1437,6 +1414,8 @@ int main(int argc, char *argv[])
|
||||
// Calculate a normal per region
|
||||
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
vectorField regionNormals(regionPoints.size(), vector::zero);
|
||||
vectorField regionCentres(regionPoints.size(), vector::zero);
|
||||
labelList nRegionFaces(regionPoints.size(), 0);
|
||||
|
||||
forAll(pointRegions, faceI)
|
||||
{
|
||||
@ -1446,9 +1425,15 @@ int main(int argc, char *argv[])
|
||||
{
|
||||
label region = pointRegions[faceI][fp];
|
||||
regionNormals[region] += extrudePatch.faceNormals()[faceI];
|
||||
regionCentres[region] += extrudePatch.faceCentres()[faceI];
|
||||
nRegionFaces[region]++;
|
||||
}
|
||||
}
|
||||
regionNormals /= mag(regionNormals);
|
||||
forAll(regionCentres, regionI)
|
||||
{
|
||||
regionCentres[regionI] /= nRegionFaces[regionI];
|
||||
}
|
||||
|
||||
|
||||
// Constrain&sync normals on points that are on coupled patches.
|
||||
@ -1463,10 +1448,13 @@ int main(int argc, char *argv[])
|
||||
);
|
||||
|
||||
// For debugging: dump hedgehog plot of normals
|
||||
if (false)
|
||||
{
|
||||
OFstream str(runTime.path()/"regionNormals.obj");
|
||||
label vertI = 0;
|
||||
|
||||
scalar thickness = model().sumThickness(1);
|
||||
|
||||
forAll(pointRegions, faceI)
|
||||
{
|
||||
const face& f = extrudeFaces[faceI];
|
||||
@ -1486,6 +1474,16 @@ int main(int argc, char *argv[])
|
||||
}
|
||||
|
||||
|
||||
// Use model to create displacements of first layer
|
||||
vectorField firstDisp(regionNormals.size());
|
||||
forAll(firstDisp, regionI)
|
||||
{
|
||||
const point& regionPt = regionCentres[regionI];
|
||||
const vector& n = regionNormals[regionI];
|
||||
firstDisp[regionI] = model()(regionPt, n, 1) - regionPt;
|
||||
}
|
||||
|
||||
|
||||
|
||||
// Create a new mesh
|
||||
// ~~~~~~~~~~~~~~~~~
|
||||
@ -1500,7 +1498,9 @@ int main(int argc, char *argv[])
|
||||
|
||||
extruder.setRefinement
|
||||
(
|
||||
thickness*regionNormals,
|
||||
firstDisp, // first displacement
|
||||
model().expansionRatio(),
|
||||
model().nLayers(), // nLayers
|
||||
extrudeTopPatchID,
|
||||
extrudeBottomPatchID,
|
||||
extrudeEdgePatches,
|
||||
@ -0,0 +1,78 @@
|
||||
/*--------------------------------*- C++ -*----------------------------------*\
|
||||
| ========= | |
|
||||
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
|
||||
| \\ / O peration | Version: dev |
|
||||
| \\ / A nd | Web: www.OpenFOAM.com |
|
||||
| \\/ M anipulation | |
|
||||
\*---------------------------------------------------------------------------*/
|
||||
FoamFile
|
||||
{
|
||||
version 2.0;
|
||||
format ascii;
|
||||
class dictionary;
|
||||
object extrudeToRegionMeshDict;
|
||||
}
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
// Name of region to create
|
||||
region liquidFilm;
|
||||
|
||||
// faceZones to extrude
|
||||
faceZones (f0);
|
||||
|
||||
// Extrude 1D-columns of cells?
|
||||
oneD false;
|
||||
|
||||
//- Extrusion model to use. The only logical choice is linearNormal?
|
||||
|
||||
//- Linear extrusion in normal direction
|
||||
extrudeModel linearNormal;
|
||||
|
||||
//- Linear extrusion in specified direction
|
||||
//extrudeModel linearDirection;
|
||||
|
||||
//- Wedge extrusion. If nLayers is 1 assumes symmetry around plane.
|
||||
// extrudeModel wedge;
|
||||
|
||||
//- Extrudes into sphere around (0 0 0)
|
||||
//extrudeModel linearRadial;
|
||||
|
||||
//- Extrudes into sphere with grading according to pressure (atmospherics)
|
||||
//extrudeModel sigmaRadial;
|
||||
|
||||
nLayers 10;
|
||||
|
||||
expansionRatio 0.9;
|
||||
|
||||
linearNormalCoeffs
|
||||
{
|
||||
thickness 0.05;
|
||||
}
|
||||
|
||||
wedgeCoeffs
|
||||
{
|
||||
axisPt (0 0.1 -0.05);
|
||||
axis (-1 0 0);
|
||||
angle 360; // For nLayers=1 assume symmetry so angle/2 on each side
|
||||
}
|
||||
|
||||
linearDirectionCoeffs
|
||||
{
|
||||
direction (0 1 0);
|
||||
thickness 0.05;
|
||||
}
|
||||
|
||||
linearRadialCoeffs
|
||||
{
|
||||
R 0.1;
|
||||
}
|
||||
|
||||
sigmaRadialCoeffs
|
||||
{
|
||||
RTbyg 1;
|
||||
pRef 1;
|
||||
pStrat 1;
|
||||
}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
@ -55,23 +55,6 @@ namespace Foam
|
||||
defineTemplateTypeNameAndDebug(IOPtrList<dictionary>, 0);
|
||||
}
|
||||
|
||||
// Combine operator to synchronise points. We choose point nearest to origin so
|
||||
// we can use e.g. great,great,great as null value.
|
||||
class nearestEqOp
|
||||
{
|
||||
|
||||
public:
|
||||
|
||||
void operator()(vector& x, const vector& y) const
|
||||
{
|
||||
if (magSqr(y) < magSqr(x))
|
||||
{
|
||||
x = y;
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
void changePatchID
|
||||
(
|
||||
const polyMesh& mesh,
|
||||
@ -854,7 +837,7 @@ int main(int argc, char *argv[])
|
||||
(
|
||||
mesh,
|
||||
newPoints,
|
||||
nearestEqOp(),
|
||||
minMagSqrEqOp<vector>(),
|
||||
point(GREAT, GREAT, GREAT)
|
||||
);
|
||||
|
||||
|
||||
@ -36,12 +36,12 @@ using namespace Foam;
|
||||
template<class Type>
|
||||
void ensightCloudField
|
||||
(
|
||||
const Foam::IOobject& fieldObject,
|
||||
const Foam::fileName& postProcPath,
|
||||
const Foam::word& prepend,
|
||||
const Foam::label timeIndex,
|
||||
const Foam::word& cloudName,
|
||||
Foam::Ostream& ensightCaseFile,
|
||||
const IOobject& fieldObject,
|
||||
const fileName& postProcPath,
|
||||
const word& prepend,
|
||||
const label timeIndex,
|
||||
const word& cloudName,
|
||||
Ostream& ensightCaseFile,
|
||||
const bool dataExists
|
||||
)
|
||||
{
|
||||
|
||||
@ -42,13 +42,13 @@ SourceFiles
|
||||
template<class Type>
|
||||
void ensightCloudField
|
||||
(
|
||||
const Foam::IOobject& fieldObject,
|
||||
const Foam::fileName& postProcPath,
|
||||
const Foam::word& prepend,
|
||||
const Foam::label timeIndex,
|
||||
const Foam::word& timeFile,
|
||||
const Foam::word& cloudName,
|
||||
Foam::Ostream& ensightCaseFile,
|
||||
const IOobject& fieldObject,
|
||||
const fileName& postProcPath,
|
||||
const word& prepend,
|
||||
const label timeIndex,
|
||||
const word& timeFile,
|
||||
const word& cloudName,
|
||||
Ostream& ensightCaseFile,
|
||||
const bool dataExists
|
||||
);
|
||||
|
||||
|
||||
@ -105,11 +105,11 @@ void writeField
|
||||
template<class Type>
|
||||
bool writePatchField
|
||||
(
|
||||
const Foam::Field<Type>& pf,
|
||||
const Foam::label patchi,
|
||||
const Foam::label ensightPatchI,
|
||||
const Foam::faceSets& boundaryFaceSet,
|
||||
const Foam::ensightMesh::nFacePrimitives& nfp,
|
||||
const Field<Type>& pf,
|
||||
const label patchi,
|
||||
const label ensightPatchI,
|
||||
const faceSets& boundaryFaceSet,
|
||||
const ensightMesh::nFacePrimitives& nfp,
|
||||
ensightStream& ensightFile
|
||||
)
|
||||
{
|
||||
@ -153,15 +153,15 @@ bool writePatchField
|
||||
template<class Type>
|
||||
void writePatchField
|
||||
(
|
||||
const Foam::word& fieldName,
|
||||
const Foam::Field<Type>& pf,
|
||||
const Foam::word& patchName,
|
||||
const Foam::ensightMesh& eMesh,
|
||||
const Foam::fileName& postProcPath,
|
||||
const Foam::word& prepend,
|
||||
const Foam::label timeIndex,
|
||||
const word& fieldName,
|
||||
const Field<Type>& pf,
|
||||
const word& patchName,
|
||||
const ensightMesh& eMesh,
|
||||
const fileName& postProcPath,
|
||||
const word& prepend,
|
||||
const label timeIndex,
|
||||
const bool binary,
|
||||
Foam::Ostream& ensightCaseFile
|
||||
Ostream& ensightCaseFile
|
||||
)
|
||||
{
|
||||
const Time& runTime = eMesh.mesh().time();
|
||||
@ -271,12 +271,12 @@ template<class Type>
|
||||
void ensightField
|
||||
(
|
||||
const GeometricField<Type, fvPatchField, volMesh>& vf,
|
||||
const Foam::ensightMesh& eMesh,
|
||||
const Foam::fileName& postProcPath,
|
||||
const Foam::word& prepend,
|
||||
const Foam::label timeIndex,
|
||||
const ensightMesh& eMesh,
|
||||
const fileName& postProcPath,
|
||||
const word& prepend,
|
||||
const label timeIndex,
|
||||
const bool binary,
|
||||
Foam::Ostream& ensightCaseFile
|
||||
Ostream& ensightCaseFile
|
||||
)
|
||||
{
|
||||
Info<< "Converting field " << vf.name() << endl;
|
||||
@ -500,12 +500,12 @@ template<class Type>
|
||||
void ensightPointField
|
||||
(
|
||||
const GeometricField<Type, pointPatchField, pointMesh>& pf,
|
||||
const Foam::ensightMesh& eMesh,
|
||||
const Foam::fileName& postProcPath,
|
||||
const Foam::word& prepend,
|
||||
const Foam::label timeIndex,
|
||||
const ensightMesh& eMesh,
|
||||
const fileName& postProcPath,
|
||||
const word& prepend,
|
||||
const label timeIndex,
|
||||
const bool binary,
|
||||
Foam::Ostream& ensightCaseFile
|
||||
Ostream& ensightCaseFile
|
||||
)
|
||||
{
|
||||
Info<< "Converting field " << pf.name() << endl;
|
||||
@ -679,14 +679,14 @@ void ensightPointField
|
||||
template<class Type>
|
||||
void ensightField
|
||||
(
|
||||
const Foam::IOobject& fieldObject,
|
||||
const Foam::ensightMesh& eMesh,
|
||||
const Foam::fileName& postProcPath,
|
||||
const Foam::word& prepend,
|
||||
const Foam::label timeIndex,
|
||||
const IOobject& fieldObject,
|
||||
const ensightMesh& eMesh,
|
||||
const fileName& postProcPath,
|
||||
const word& prepend,
|
||||
const label timeIndex,
|
||||
const bool binary,
|
||||
const bool nodeValues,
|
||||
Foam::Ostream& ensightCaseFile
|
||||
Ostream& ensightCaseFile
|
||||
)
|
||||
{
|
||||
// Read field
|
||||
|
||||
@ -37,10 +37,10 @@ using namespace Foam;
|
||||
|
||||
void ensightParticlePositions
|
||||
(
|
||||
const Foam::fvMesh& mesh,
|
||||
const Foam::fileName& postProcPath,
|
||||
const Foam::word& timeFile,
|
||||
const Foam::word& cloudName,
|
||||
const fvMesh& mesh,
|
||||
const fileName& postProcPath,
|
||||
const word& timeFile,
|
||||
const word& cloudName,
|
||||
const bool dataExists
|
||||
)
|
||||
{
|
||||
|
||||
@ -33,9 +33,6 @@ Description
|
||||
|
||||
#include "argList.H"
|
||||
#include "ListOps.H"
|
||||
#include "face.H"
|
||||
#include "tetPointRef.H"
|
||||
#include "triFaceList.H"
|
||||
#include "triSurface.H"
|
||||
#include "OFstream.H"
|
||||
#include "meshTools.H"
|
||||
@ -43,242 +40,12 @@ Description
|
||||
#include "transform.H"
|
||||
#include "IOmanip.H"
|
||||
#include "Pair.H"
|
||||
#include "momentOfInertia.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
using namespace Foam;
|
||||
|
||||
void massPropertiesSolid
|
||||
(
|
||||
const pointField& pts,
|
||||
const triFaceList& triFaces,
|
||||
scalar density,
|
||||
scalar& mass,
|
||||
vector& cM,
|
||||
tensor& J
|
||||
)
|
||||
{
|
||||
// Reimplemented from: Wm4PolyhedralMassProperties.cpp
|
||||
// File Version: 4.10.0 (2009/11/18)
|
||||
|
||||
// Geometric Tools, LC
|
||||
// Copyright (c) 1998-2010
|
||||
// Distributed under the Boost Software License, Version 1.0.
|
||||
// http://www.boost.org/LICENSE_1_0.txt
|
||||
// http://www.geometrictools.com/License/Boost/LICENSE_1_0.txt
|
||||
|
||||
// Boost Software License - Version 1.0 - August 17th, 2003
|
||||
|
||||
// Permission is hereby granted, free of charge, to any person or
|
||||
// organization obtaining a copy of the software and accompanying
|
||||
// documentation covered by this license (the "Software") to use,
|
||||
// reproduce, display, distribute, execute, and transmit the
|
||||
// Software, and to prepare derivative works of the Software, and
|
||||
// to permit third-parties to whom the Software is furnished to do
|
||||
// so, all subject to the following:
|
||||
|
||||
// The copyright notices in the Software and this entire
|
||||
// statement, including the above license grant, this restriction
|
||||
// and the following disclaimer, must be included in all copies of
|
||||
// the Software, in whole or in part, and all derivative works of
|
||||
// the Software, unless such copies or derivative works are solely
|
||||
// in the form of machine-executable object code generated by a
|
||||
// source language processor.
|
||||
|
||||
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
|
||||
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
|
||||
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE AND
|
||||
// NON-INFRINGEMENT. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR
|
||||
// ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE FOR ANY DAMAGES OR
|
||||
// OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
|
||||
// ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE
|
||||
// USE OR OTHER DEALINGS IN THE SOFTWARE.
|
||||
|
||||
const scalar r6 = 1.0/6.0;
|
||||
const scalar r24 = 1.0/24.0;
|
||||
const scalar r60 = 1.0/60.0;
|
||||
const scalar r120 = 1.0/120.0;
|
||||
|
||||
// order: 1, x, y, z, x^2, y^2, z^2, xy, yz, zx
|
||||
scalarField integrals(10, 0.0);
|
||||
|
||||
forAll(triFaces, i)
|
||||
{
|
||||
const triFace& tri(triFaces[i]);
|
||||
|
||||
// vertices of triangle i
|
||||
vector v0 = pts[tri[0]];
|
||||
vector v1 = pts[tri[1]];
|
||||
vector v2 = pts[tri[2]];
|
||||
|
||||
// cross product of edges
|
||||
vector eA = v1 - v0;
|
||||
vector eB = v2 - v0;
|
||||
vector n = eA ^ eB;
|
||||
|
||||
// compute integral terms
|
||||
scalar tmp0, tmp1, tmp2;
|
||||
|
||||
scalar f1x, f2x, f3x, g0x, g1x, g2x;
|
||||
|
||||
tmp0 = v0.x() + v1.x();
|
||||
f1x = tmp0 + v2.x();
|
||||
tmp1 = v0.x()*v0.x();
|
||||
tmp2 = tmp1 + v1.x()*tmp0;
|
||||
f2x = tmp2 + v2.x()*f1x;
|
||||
f3x = v0.x()*tmp1 + v1.x()*tmp2 + v2.x()*f2x;
|
||||
g0x = f2x + v0.x()*(f1x + v0.x());
|
||||
g1x = f2x + v1.x()*(f1x + v1.x());
|
||||
g2x = f2x + v2.x()*(f1x + v2.x());
|
||||
|
||||
scalar f1y, f2y, f3y, g0y, g1y, g2y;
|
||||
|
||||
tmp0 = v0.y() + v1.y();
|
||||
f1y = tmp0 + v2.y();
|
||||
tmp1 = v0.y()*v0.y();
|
||||
tmp2 = tmp1 + v1.y()*tmp0;
|
||||
f2y = tmp2 + v2.y()*f1y;
|
||||
f3y = v0.y()*tmp1 + v1.y()*tmp2 + v2.y()*f2y;
|
||||
g0y = f2y + v0.y()*(f1y + v0.y());
|
||||
g1y = f2y + v1.y()*(f1y + v1.y());
|
||||
g2y = f2y + v2.y()*(f1y + v2.y());
|
||||
|
||||
scalar f1z, f2z, f3z, g0z, g1z, g2z;
|
||||
|
||||
tmp0 = v0.z() + v1.z();
|
||||
f1z = tmp0 + v2.z();
|
||||
tmp1 = v0.z()*v0.z();
|
||||
tmp2 = tmp1 + v1.z()*tmp0;
|
||||
f2z = tmp2 + v2.z()*f1z;
|
||||
f3z = v0.z()*tmp1 + v1.z()*tmp2 + v2.z()*f2z;
|
||||
g0z = f2z + v0.z()*(f1z + v0.z());
|
||||
g1z = f2z + v1.z()*(f1z + v1.z());
|
||||
g2z = f2z + v2.z()*(f1z + v2.z());
|
||||
|
||||
// update integrals
|
||||
integrals[0] += n.x()*f1x;
|
||||
integrals[1] += n.x()*f2x;
|
||||
integrals[2] += n.y()*f2y;
|
||||
integrals[3] += n.z()*f2z;
|
||||
integrals[4] += n.x()*f3x;
|
||||
integrals[5] += n.y()*f3y;
|
||||
integrals[6] += n.z()*f3z;
|
||||
integrals[7] += n.x()*(v0.y()*g0x + v1.y()*g1x + v2.y()*g2x);
|
||||
integrals[8] += n.y()*(v0.z()*g0y + v1.z()*g1y + v2.z()*g2y);
|
||||
integrals[9] += n.z()*(v0.x()*g0z + v1.x()*g1z + v2.x()*g2z);
|
||||
}
|
||||
|
||||
integrals[0] *= r6;
|
||||
integrals[1] *= r24;
|
||||
integrals[2] *= r24;
|
||||
integrals[3] *= r24;
|
||||
integrals[4] *= r60;
|
||||
integrals[5] *= r60;
|
||||
integrals[6] *= r60;
|
||||
integrals[7] *= r120;
|
||||
integrals[8] *= r120;
|
||||
integrals[9] *= r120;
|
||||
|
||||
// mass
|
||||
mass = integrals[0];
|
||||
|
||||
// center of mass
|
||||
cM = vector(integrals[1], integrals[2], integrals[3])/mass;
|
||||
|
||||
// inertia relative to origin
|
||||
J.xx() = integrals[5] + integrals[6];
|
||||
J.xy() = -integrals[7];
|
||||
J.xz() = -integrals[9];
|
||||
J.yx() = J.xy();
|
||||
J.yy() = integrals[4] + integrals[6];
|
||||
J.yz() = -integrals[8];
|
||||
J.zx() = J.xz();
|
||||
J.zy() = J.yz();
|
||||
J.zz() = integrals[4] + integrals[5];
|
||||
|
||||
// inertia relative to center of mass
|
||||
J -= mass*((cM & cM)*I - cM*cM);
|
||||
|
||||
// Apply density
|
||||
mass *= density;
|
||||
J *= density;
|
||||
}
|
||||
|
||||
|
||||
void massPropertiesShell
|
||||
(
|
||||
const pointField& pts,
|
||||
const triFaceList& triFaces,
|
||||
scalar density,
|
||||
scalar& mass,
|
||||
vector& cM,
|
||||
tensor& J
|
||||
)
|
||||
{
|
||||
// Reset properties for accumulation
|
||||
|
||||
mass = 0.0;
|
||||
cM = vector::zero;
|
||||
J = tensor::zero;
|
||||
|
||||
// Find centre of mass
|
||||
|
||||
forAll(triFaces, i)
|
||||
{
|
||||
const triFace& tri(triFaces[i]);
|
||||
|
||||
triPointRef t
|
||||
(
|
||||
pts[tri[0]],
|
||||
pts[tri[1]],
|
||||
pts[tri[2]]
|
||||
);
|
||||
|
||||
scalar triMag = t.mag();
|
||||
|
||||
cM += triMag*t.centre();
|
||||
|
||||
mass += triMag;
|
||||
}
|
||||
|
||||
cM /= mass;
|
||||
|
||||
mass *= density;
|
||||
|
||||
// Find inertia around centre of mass
|
||||
|
||||
forAll(triFaces, i)
|
||||
{
|
||||
const triFace& tri(triFaces[i]);
|
||||
|
||||
J += triPointRef
|
||||
(
|
||||
pts[tri[0]],
|
||||
pts[tri[1]],
|
||||
pts[tri[2]]
|
||||
).inertia(cM, density);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
tensor applyParallelAxisTheorem
|
||||
(
|
||||
scalar m,
|
||||
const vector& cM,
|
||||
const tensor& J,
|
||||
const vector& refPt
|
||||
)
|
||||
{
|
||||
// The displacement vector (refPt = cM) is the displacement of the
|
||||
// new reference point from the centre of mass of the body that
|
||||
// the inertia tensor applies to.
|
||||
|
||||
vector d = (refPt - cM);
|
||||
|
||||
return J + m*((d & d)*I - d*d);
|
||||
}
|
||||
|
||||
|
||||
int main(int argc, char *argv[])
|
||||
{
|
||||
argList::addNote
|
||||
@ -321,40 +88,17 @@ int main(int argc, char *argv[])
|
||||
|
||||
triSurface surf(surfFileName);
|
||||
|
||||
triFaceList faces(surf.size());
|
||||
|
||||
forAll(surf, i)
|
||||
{
|
||||
faces[i] = triFace(surf[i]);
|
||||
}
|
||||
|
||||
scalar m = 0.0;
|
||||
vector cM = vector::zero;
|
||||
tensor J = tensor::zero;
|
||||
|
||||
if (args.optionFound("shellProperties"))
|
||||
{
|
||||
massPropertiesShell
|
||||
(
|
||||
surf.points(),
|
||||
faces,
|
||||
density,
|
||||
m,
|
||||
cM,
|
||||
J
|
||||
);
|
||||
momentOfInertia::massPropertiesShell(surf, density, m, cM, J);
|
||||
}
|
||||
else
|
||||
{
|
||||
massPropertiesSolid
|
||||
(
|
||||
surf.points(),
|
||||
faces,
|
||||
density,
|
||||
m,
|
||||
cM,
|
||||
J
|
||||
);
|
||||
momentOfInertia::massPropertiesSolid(surf, density, m, cM, J);
|
||||
}
|
||||
|
||||
if (m < 0)
|
||||
@ -583,7 +327,7 @@ int main(int argc, char *argv[])
|
||||
showTransform = false;
|
||||
}
|
||||
|
||||
Info<< nl << setprecision(10)
|
||||
Info<< nl << setprecision(12)
|
||||
<< "Density: " << density << nl
|
||||
<< "Mass: " << m << nl
|
||||
<< "Centre of mass: " << cM << nl
|
||||
@ -615,7 +359,7 @@ int main(int argc, char *argv[])
|
||||
if (calcAroundRefPt)
|
||||
{
|
||||
Info<< nl << "Inertia tensor relative to " << refPt << ": " << nl
|
||||
<< applyParallelAxisTheorem(m, cM, J, refPt)
|
||||
<< momentOfInertia::applyParallelAxisTheorem(m, cM, J, refPt)
|
||||
<< endl;
|
||||
}
|
||||
|
||||
|
||||
@ -41,6 +41,7 @@ SourceFiles
|
||||
#include "tetPointRef.H"
|
||||
#include "triPointRef.H"
|
||||
#include "polyMesh.H"
|
||||
#include "triFace.H"
|
||||
#include "face.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
@ -146,6 +147,10 @@ public:
|
||||
// mesh face for this tet from the supplied mesh
|
||||
inline triPointRef faceTri(const polyMesh& mesh) const;
|
||||
|
||||
//- Return the point indices corresponding to the tri on the mesh
|
||||
// face for this tet from the supplied mesh
|
||||
inline triFace faceTriIs(const polyMesh& mesh) const;
|
||||
|
||||
//- Return the geometry corresponding to the tri on the
|
||||
// mesh face for this tet from the supplied mesh using
|
||||
// the old position
|
||||
|
||||
@ -122,6 +122,21 @@ Foam::triPointRef Foam::tetIndices::faceTri(const polyMesh& mesh) const
|
||||
}
|
||||
|
||||
|
||||
Foam::triFace Foam::tetIndices::faceTriIs(const polyMesh& mesh) const
|
||||
{
|
||||
const faceList& pFaces = mesh.faces();
|
||||
|
||||
const Foam::face& f = pFaces[faceI_];
|
||||
|
||||
return triFace
|
||||
(
|
||||
f[faceBasePtI_],
|
||||
f[facePtAI_],
|
||||
f[facePtBI_]
|
||||
);
|
||||
}
|
||||
|
||||
|
||||
Foam::triPointRef Foam::tetIndices::oldFaceTri(const polyMesh& mesh) const
|
||||
{
|
||||
const pointField& oldPPts = mesh.oldPoints();
|
||||
|
||||
@ -156,29 +156,12 @@ inline Point tetrahedron<Point, PointRef>::circumCentre() const
|
||||
|
||||
if (Foam::mag(denom) < ROOTVSMALL)
|
||||
{
|
||||
// Degenerate tet. Use test of the individual triangles.
|
||||
{
|
||||
point triCentre = triPointRef(a_, b_, c_).circumCentre();
|
||||
if (magSqr(d_ - triCentre) < magSqr(a_ - triCentre))
|
||||
{
|
||||
return triCentre;
|
||||
}
|
||||
}
|
||||
{
|
||||
point triCentre = triPointRef(a_, b_, d_).circumCentre();
|
||||
if (magSqr(c_ - triCentre) < magSqr(a_ - triCentre))
|
||||
{
|
||||
return triCentre;
|
||||
}
|
||||
}
|
||||
{
|
||||
point triCentre = triPointRef(a_, c_, d_).circumCentre();
|
||||
if (magSqr(b_ - triCentre) < magSqr(a_ - triCentre))
|
||||
{
|
||||
return triCentre;
|
||||
}
|
||||
}
|
||||
return triPointRef(b_, c_, d_).circumCentre();
|
||||
WarningIn("Point tetrahedron<Point, PointRef>::circumCentre() const")
|
||||
<< "Degenerate tetrahedron:" << nl << *this << nl
|
||||
<<"Returning centre instead of circumCentre."
|
||||
<< endl;
|
||||
|
||||
return centre();
|
||||
}
|
||||
|
||||
return a_ + 0.5*(a + num/denom);
|
||||
@ -188,16 +171,43 @@ inline Point tetrahedron<Point, PointRef>::circumCentre() const
|
||||
template<class Point, class PointRef>
|
||||
inline scalar tetrahedron<Point, PointRef>::circumRadius() const
|
||||
{
|
||||
return Foam::mag(a_ - circumCentre());
|
||||
vector a = b_ - a_;
|
||||
vector b = c_ - a_;
|
||||
vector c = d_ - a_;
|
||||
|
||||
scalar lambda = magSqr(c) - (a & c);
|
||||
scalar mu = magSqr(b) - (a & b);
|
||||
|
||||
vector ba = b ^ a;
|
||||
vector ca = c ^ a;
|
||||
|
||||
vector num = lambda*ba - mu*ca;
|
||||
scalar denom = (c & ba);
|
||||
|
||||
if (Foam::mag(denom) < ROOTVSMALL)
|
||||
{
|
||||
WarningIn("Point tetrahedron<Point, PointRef>::circumCentre() const")
|
||||
<< "Degenerate tetrahedron:" << nl << *this << nl
|
||||
<< "Returning GREAT for circumRadius."
|
||||
<< endl;
|
||||
|
||||
return GREAT;
|
||||
}
|
||||
|
||||
return Foam::mag(0.5*(a + num/denom));
|
||||
}
|
||||
|
||||
|
||||
template<class Point, class PointRef>
|
||||
inline scalar tetrahedron<Point, PointRef>::quality() const
|
||||
{
|
||||
// Note: 8/(9*sqrt(3)) = 0.5132002393
|
||||
|
||||
return mag()/(0.5132002393*pow3(min(circumRadius(), GREAT)) + ROOTVSMALL);
|
||||
return
|
||||
mag()
|
||||
/(
|
||||
8.0/(9.0*sqrt(3.0))
|
||||
*pow3(min(circumRadius(), GREAT))
|
||||
+ ROOTVSMALL
|
||||
);
|
||||
}
|
||||
|
||||
|
||||
@ -266,7 +276,8 @@ scalar tetrahedron<Point, PointRef>::barycentric
|
||||
"const point& pt"
|
||||
") const"
|
||||
)
|
||||
<< "Degenerate tetrahedron - returning 1/4 barycentric coordinates."
|
||||
<< "Degenerate tetrahedron:" << nl << *this << nl
|
||||
<< "Returning 1/4 barycentric coordinates."
|
||||
<< endl;
|
||||
|
||||
bary = List<scalar>(4, 0.25);
|
||||
|
||||
@ -109,9 +109,9 @@ inline vector triangle<Point, PointRef>::normal() const
|
||||
template<class Point, class PointRef>
|
||||
inline Point triangle<Point, PointRef>::circumCentre() const
|
||||
{
|
||||
scalar d1 = (c_ - a_)&(b_ - a_);
|
||||
scalar d2 = -(c_ - b_)&(b_ - a_);
|
||||
scalar d3 = (c_ - a_)&(c_ - b_);
|
||||
scalar d1 = (c_ - a_) & (b_ - a_);
|
||||
scalar d2 = -(c_ - b_) & (b_ - a_);
|
||||
scalar d3 = (c_ - a_) & (c_ - b_);
|
||||
|
||||
scalar c1 = d2*d3;
|
||||
scalar c2 = d3*d1;
|
||||
@ -119,6 +119,16 @@ inline Point triangle<Point, PointRef>::circumCentre() const
|
||||
|
||||
scalar c = c1 + c2 + c3;
|
||||
|
||||
if (Foam::mag(c) < ROOTVSMALL)
|
||||
{
|
||||
WarningIn("Point triangle<Point, PointRef>::circumCentre() const")
|
||||
<< "Degenerate triangle:" << nl << *this << nl
|
||||
<< "Returning centre instead of circumCentre."
|
||||
<< endl;
|
||||
|
||||
return centre();
|
||||
}
|
||||
|
||||
return
|
||||
(
|
||||
((c2 + c3)*a_ + (c3 + c1)*b_ + (c1 + c2)*c_)/(2*c)
|
||||
@ -129,14 +139,19 @@ inline Point triangle<Point, PointRef>::circumCentre() const
|
||||
template<class Point, class PointRef>
|
||||
inline scalar triangle<Point, PointRef>::circumRadius() const
|
||||
{
|
||||
scalar d1 = (c_ - a_) & (b_ - a_);
|
||||
scalar d2 = - (c_ - b_) & (b_ - a_);
|
||||
scalar d3 = (c_ - a_) & (c_ - b_);
|
||||
scalar d1 = (c_ - a_) & (b_ - a_);
|
||||
scalar d2 = -(c_ - b_) & (b_ - a_);
|
||||
scalar d3 = (c_ - a_) & (c_ - b_);
|
||||
|
||||
scalar denom = d2*d3 + d3*d1 + d1*d2;
|
||||
|
||||
if (Foam::mag(denom) < VSMALL)
|
||||
{
|
||||
WarningIn("scalar triangle<Point, PointRef>::circumRadius() const")
|
||||
<< "Degenerate triangle:" << nl << *this << nl
|
||||
<< "Returning GREAT for circumRadius."
|
||||
<< endl;
|
||||
|
||||
return GREAT;
|
||||
}
|
||||
else
|
||||
@ -151,16 +166,7 @@ inline scalar triangle<Point, PointRef>::circumRadius() const
|
||||
template<class Point, class PointRef>
|
||||
inline scalar triangle<Point, PointRef>::quality() const
|
||||
{
|
||||
// Note: 3*sqr(3)/(4*pi) = 0.4134966716
|
||||
|
||||
return
|
||||
mag()
|
||||
/ (
|
||||
constant::mathematical::pi
|
||||
*Foam::sqr(circumRadius())
|
||||
*0.4134966716
|
||||
+ VSMALL
|
||||
);
|
||||
return mag()/(Foam::sqr(circumRadius())*3.0*sqrt(3.0)/4.0 + VSMALL);
|
||||
}
|
||||
|
||||
|
||||
@ -264,7 +270,8 @@ scalar triangle<Point, PointRef>::barycentric
|
||||
"const point& pt"
|
||||
") const"
|
||||
)
|
||||
<< "Degenerate triangle - returning 1/3 barycentric coordinates."
|
||||
<< "Degenerate triangle:" << nl << *this << nl
|
||||
<< "Returning 1/3 barycentric coordinates."
|
||||
<< endl;
|
||||
|
||||
bary = List<scalar>(3, 1.0/3.0);
|
||||
@ -490,7 +497,7 @@ pointHit triangle<Point, PointRef>::nearestPointClassify
|
||||
) const
|
||||
{
|
||||
// Adapted from:
|
||||
// Real-time collision detection, Christer Ericson, 2005, 136-142
|
||||
// Real-time collision detection, Christer Ericson, 2005, p136-142
|
||||
|
||||
// Check if P in vertex region outside A
|
||||
vector ab = b_ - a_;
|
||||
@ -528,6 +535,27 @@ pointHit triangle<Point, PointRef>::nearestPointClassify
|
||||
|
||||
if (vc <= 0.0 && d1 >= 0.0 && d3 <= 0.0)
|
||||
{
|
||||
if ((d1 - d3) < ROOTVSMALL)
|
||||
{
|
||||
WarningIn
|
||||
(
|
||||
"pointHit triangle<Point, PointRef>::nearestPointClassify"
|
||||
"("
|
||||
"const point& p,"
|
||||
"label& nearType,"
|
||||
"label& nearLabel"
|
||||
") const"
|
||||
)
|
||||
<< "Degenerate triangle:" << nl << *this << nl
|
||||
<< "d1, d3: " << d1 << ", " << d3 << endl;
|
||||
|
||||
// For d1 = d3, a_ and b_ are likely coincident.
|
||||
|
||||
nearType = POINT;
|
||||
nearLabel = 0;
|
||||
return pointHit(false, a_, Foam::mag(a_ - p), true);
|
||||
}
|
||||
|
||||
// barycentric coordinates (1-v, v, 0)
|
||||
scalar v = d1/(d1 - d3);
|
||||
|
||||
@ -556,6 +584,27 @@ pointHit triangle<Point, PointRef>::nearestPointClassify
|
||||
|
||||
if (vb <= 0.0 && d2 >= 0.0 && d6 <= 0.0)
|
||||
{
|
||||
if ((d2 - d6) < ROOTVSMALL)
|
||||
{
|
||||
WarningIn
|
||||
(
|
||||
"pointHit triangle<Point, PointRef>::nearestPointClassify"
|
||||
"("
|
||||
"const point& p,"
|
||||
"label& nearType,"
|
||||
"label& nearLabel"
|
||||
") const"
|
||||
)
|
||||
<< "Degenerate triangle:" << nl << *this << nl
|
||||
<< "d2, d6: " << d2 << ", " << d6 << endl;
|
||||
|
||||
// For d2 = d6, a_ and c_ are likely coincident.
|
||||
|
||||
nearType = POINT;
|
||||
nearLabel = 0;
|
||||
return pointHit(false, a_, Foam::mag(a_ - p), true);
|
||||
}
|
||||
|
||||
// barycentric coordinates (1-w, 0, w)
|
||||
scalar w = d2/(d2 - d6);
|
||||
|
||||
@ -570,6 +619,28 @@ pointHit triangle<Point, PointRef>::nearestPointClassify
|
||||
|
||||
if (va <= 0.0 && (d4 - d3) >= 0.0 && (d5 - d6) >= 0.0)
|
||||
{
|
||||
if (((d4 - d3) + (d5 - d6)) < ROOTVSMALL)
|
||||
{
|
||||
WarningIn
|
||||
(
|
||||
"pointHit triangle<Point, PointRef>::nearestPointClassify"
|
||||
"("
|
||||
"const point& p,"
|
||||
"label& nearType,"
|
||||
"label& nearLabel"
|
||||
") const"
|
||||
)
|
||||
<< "Degenerate triangle:" << nl << *this << nl
|
||||
<< "(d4 - d3), (d6 - d5): " << (d4 - d3) << ", " << (d6 - d5)
|
||||
<< endl;
|
||||
|
||||
// For (d4 - d3) = (d6 - d5), b_ and c_ are likely coincident.
|
||||
|
||||
nearType = POINT;
|
||||
nearLabel = 1;
|
||||
return pointHit(false, b_, Foam::mag(b_ - p), true);
|
||||
}
|
||||
|
||||
// barycentric coordinates (0, 1-w, w)
|
||||
scalar w = (d4 - d3)/((d4 - d3) + (d5 - d6));
|
||||
|
||||
@ -581,6 +652,28 @@ pointHit triangle<Point, PointRef>::nearestPointClassify
|
||||
|
||||
// P inside face region. Compute Q through its barycentric
|
||||
// coordinates (u, v, w)
|
||||
|
||||
if ((va + vb + vc) < ROOTVSMALL)
|
||||
{
|
||||
WarningIn
|
||||
(
|
||||
"pointHit triangle<Point, PointRef>::nearestPointClassify"
|
||||
"("
|
||||
"const point& p,"
|
||||
"label& nearType,"
|
||||
"label& nearLabel"
|
||||
") const"
|
||||
)
|
||||
<< "Degenerate triangle:" << nl << *this << nl
|
||||
<< "va, vb, vc: " << va << ", " << vb << ", " << vc
|
||||
<< endl;
|
||||
|
||||
point nearPt = centre();
|
||||
nearType = NONE,
|
||||
nearLabel = -1;
|
||||
return pointHit(true, nearPt, Foam::mag(nearPt - p), false);
|
||||
}
|
||||
|
||||
scalar denom = 1.0/(va + vb + vc);
|
||||
scalar v = vb * denom;
|
||||
scalar w = vc * denom;
|
||||
|
||||
@ -76,6 +76,8 @@ EqOp(eqMag, x = mag(y))
|
||||
EqOp(plusEqMagSqr, x += magSqr(y))
|
||||
EqOp(maxEq, x = max(x, y))
|
||||
EqOp(minEq, x = min(x, y))
|
||||
EqOp(minMagSqrEq, x = (magSqr(x)<=magSqr(y) ? x : y))
|
||||
EqOp(maxMagSqrEq, x = (magSqr(x)>=magSqr(y) ? x : y))
|
||||
EqOp(andEq, x = (x && y))
|
||||
EqOp(orEq, x = (x || y))
|
||||
|
||||
|
||||
@ -48,9 +48,11 @@ inline tensor rotationTensor
|
||||
const vector& n2
|
||||
)
|
||||
{
|
||||
const scalar s = n1 & n2;
|
||||
const vector n3 = n1 ^ n2;
|
||||
return
|
||||
(n1 & n2)*I
|
||||
+ (1 - (n1 & n2))*sqr(n1 ^ n2)/(magSqr(n1 ^ n2) + VSMALL)
|
||||
s*I
|
||||
+ (1 - s)*sqr(n3)/(magSqr(n3) + VSMALL)
|
||||
+ (n2*n1 - n1*n2);
|
||||
}
|
||||
|
||||
|
||||
@ -48,7 +48,7 @@ namespace Foam
|
||||
{
|
||||
|
||||
/*---------------------------------------------------------------------------*\
|
||||
Class selfContainedDirectMappedFixedValueFvPatch Declaration
|
||||
Class selfContainedDirectMappedFixedValueFvPatchField Declaration
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
template<class Type>
|
||||
|
||||
@ -127,6 +127,7 @@ $(cellZoneSources)/setToCellZone/setToCellZone.C
|
||||
pointZoneSources = sets/pointZoneSources
|
||||
$(pointZoneSources)/setToPointZone/setToPointZone.C
|
||||
|
||||
momentOfInertia/momentOfInertia.C
|
||||
|
||||
surfaceSets/surfaceSets.C
|
||||
|
||||
|
||||
@ -21,382 +21,328 @@ License
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
Class
|
||||
momentOfInertia
|
||||
|
||||
Description
|
||||
Reimplementation of volInt.c by Brian Mirtich.
|
||||
* mirtich@cs.berkeley.edu *
|
||||
* http://www.cs.berkeley.edu/~mirtich *
|
||||
|
||||
-------------------------------------------------------------------------------
|
||||
*/
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#include "momentOfInertia.H"
|
||||
//#include "pyramidPointFaceRef.H"
|
||||
|
||||
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
|
||||
// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
|
||||
|
||||
//Foam::tensor Foam::momentOfInertia
|
||||
//(
|
||||
// const pointField& points,
|
||||
// const faceList& faces,
|
||||
// const cell& cFaces,
|
||||
// const point& cc
|
||||
//)
|
||||
//{
|
||||
// tensor t(tensor::zero);
|
||||
//
|
||||
// forAll(cFaces, i)
|
||||
// {
|
||||
// const face& f = faces[cFaces[i]];
|
||||
//
|
||||
// scalar pyrVol = pyramidPointFaceRef(f, cc).mag(points);
|
||||
//
|
||||
// vector pyrCentre = pyramidPointFaceRef(f, cc).centre(points);
|
||||
//
|
||||
// vector d = pyrCentre - cc;
|
||||
//
|
||||
// t.xx() += pyrVol*(sqr(d.y()) + sqr(d.z()));
|
||||
// t.yy() += pyrVol*(sqr(d.x()) + sqr(d.z()));
|
||||
// t.zz() += pyrVol*(sqr(d.x()) + sqr(d.y()));
|
||||
//
|
||||
// t.xy() -= pyrVol*d.x()*d.y();
|
||||
// t.xz() -= pyrVol*d.x()*d.z();
|
||||
// t.yz() -= pyrVol*d.y()*d.z();
|
||||
// }
|
||||
//
|
||||
// // Symmetric
|
||||
// t.yx() = t.xy();
|
||||
// t.zx() = t.xz();
|
||||
// t.zy() = t.yz();
|
||||
//
|
||||
// return t;
|
||||
//}
|
||||
|
||||
|
||||
#define sqr(x) ((x)*(x))
|
||||
#define pow3(x) ((x)*(x)*(x))
|
||||
|
||||
// compute various integrations over projection of face
|
||||
void Foam::compProjectionIntegrals
|
||||
void Foam::momentOfInertia::massPropertiesSolid
|
||||
(
|
||||
const pointField& points,
|
||||
const face& f,
|
||||
const direction A,
|
||||
const direction B,
|
||||
|
||||
scalar& P1,
|
||||
scalar& Pa,
|
||||
scalar& Pb,
|
||||
scalar& Paa,
|
||||
scalar& Pab,
|
||||
scalar& Pbb,
|
||||
scalar& Paaa,
|
||||
scalar& Paab,
|
||||
scalar& Pabb,
|
||||
scalar& Pbbb
|
||||
)
|
||||
{
|
||||
P1 = Pa = Pb = Paa = Pab = Pbb = Paaa = Paab = Pabb = Pbbb = 0.0;
|
||||
|
||||
forAll(f, i)
|
||||
{
|
||||
scalar a0 = points[f[i]][A];
|
||||
scalar b0 = points[f[i]][B];
|
||||
scalar a1 = points[f[(i+1) % f.size()]][A];
|
||||
scalar b1 = points[f[(i+1) % f.size()]][B];
|
||||
scalar da = a1 - a0;
|
||||
scalar db = b1 - b0;
|
||||
|
||||
scalar a0_2 = a0 * a0;
|
||||
scalar a0_3 = a0_2 * a0;
|
||||
scalar a0_4 = a0_3 * a0;
|
||||
|
||||
scalar b0_2 = b0 * b0;
|
||||
scalar b0_3 = b0_2 * b0;
|
||||
scalar b0_4 = b0_3 * b0;
|
||||
|
||||
scalar a1_2 = a1 * a1;
|
||||
scalar a1_3 = a1_2 * a1;
|
||||
|
||||
scalar b1_2 = b1 * b1;
|
||||
scalar b1_3 = b1_2 * b1;
|
||||
|
||||
scalar C1 = a1 + a0;
|
||||
|
||||
scalar Ca = a1*C1 + a0_2;
|
||||
scalar Caa = a1*Ca + a0_3;
|
||||
scalar Caaa = a1*Caa + a0_4;
|
||||
|
||||
scalar Cb = b1*(b1 + b0) + b0_2;
|
||||
scalar Cbb = b1*Cb + b0_3;
|
||||
scalar Cbbb = b1*Cbb + b0_4;
|
||||
|
||||
scalar Cab = 3*a1_2 + 2*a1*a0 + a0_2;
|
||||
scalar Kab = a1_2 + 2*a1*a0 + 3*a0_2;
|
||||
|
||||
scalar Caab = a0*Cab + 4*a1_3;
|
||||
scalar Kaab = a1*Kab + 4*a0_3;
|
||||
|
||||
scalar Cabb = 4*b1_3 + 3*b1_2*b0 + 2*b1*b0_2 + b0_3;
|
||||
scalar Kabb = b1_3 + 2*b1_2*b0 + 3*b1*b0_2 + 4*b0_3;
|
||||
|
||||
P1 += db*C1;
|
||||
Pa += db*Ca;
|
||||
Paa += db*Caa;
|
||||
Paaa += db*Caaa;
|
||||
Pb += da*Cb;
|
||||
Pbb += da*Cbb;
|
||||
Pbbb += da*Cbbb;
|
||||
Pab += db*(b1*Cab + b0*Kab);
|
||||
Paab += db*(b1*Caab + b0*Kaab);
|
||||
Pabb += da*(a1*Cabb + a0*Kabb);
|
||||
}
|
||||
|
||||
P1 /= 2.0;
|
||||
Pa /= 6.0;
|
||||
Paa /= 12.0;
|
||||
Paaa /= 20.0;
|
||||
Pb /= -6.0;
|
||||
Pbb /= -12.0;
|
||||
Pbbb /= -20.0;
|
||||
Pab /= 24.0;
|
||||
Paab /= 60.0;
|
||||
Pabb /= -60.0;
|
||||
}
|
||||
|
||||
|
||||
void Foam::compFaceIntegrals
|
||||
(
|
||||
const pointField& points,
|
||||
const face& f,
|
||||
const vector& n,
|
||||
const scalar w,
|
||||
const direction A,
|
||||
const direction B,
|
||||
const direction C,
|
||||
|
||||
scalar& Fa,
|
||||
scalar& Fb,
|
||||
scalar& Fc,
|
||||
scalar& Faa,
|
||||
scalar& Fbb,
|
||||
scalar& Fcc,
|
||||
scalar& Faaa,
|
||||
scalar& Fbbb,
|
||||
scalar& Fccc,
|
||||
scalar& Faab,
|
||||
scalar& Fbbc,
|
||||
scalar& Fcca
|
||||
)
|
||||
{
|
||||
scalar P1, Pa, Pb, Paa, Pab, Pbb, Paaa, Paab, Pabb, Pbbb;
|
||||
|
||||
compProjectionIntegrals
|
||||
(
|
||||
points,
|
||||
f,
|
||||
A,
|
||||
B,
|
||||
|
||||
P1,
|
||||
Pa,
|
||||
Pb,
|
||||
Paa,
|
||||
Pab,
|
||||
Pbb,
|
||||
Paaa,
|
||||
Paab,
|
||||
Pabb,
|
||||
Pbbb
|
||||
);
|
||||
|
||||
scalar k1 = 1 / n[C];
|
||||
scalar k2 = k1 * k1;
|
||||
scalar k3 = k2 * k1;
|
||||
scalar k4 = k3 * k1;
|
||||
|
||||
Fa = k1 * Pa;
|
||||
Fb = k1 * Pb;
|
||||
Fc = -k2 * (n[A]*Pa + n[B]*Pb + w*P1);
|
||||
|
||||
Faa = k1 * Paa;
|
||||
Fbb = k1 * Pbb;
|
||||
Fcc = k3 * (sqr(n[A])*Paa + 2*n[A]*n[B]*Pab + sqr(n[B])*Pbb
|
||||
+ w*(2*(n[A]*Pa + n[B]*Pb) + w*P1));
|
||||
|
||||
Faaa = k1 * Paaa;
|
||||
Fbbb = k1 * Pbbb;
|
||||
Fccc = -k4 * (pow3(n[A])*Paaa + 3*sqr(n[A])*n[B]*Paab
|
||||
+ 3*n[A]*sqr(n[B])*Pabb + pow3(n[B])*Pbbb
|
||||
+ 3*w*(sqr(n[A])*Paa + 2*n[A]*n[B]*Pab + sqr(n[B])*Pbb)
|
||||
+ w*w*(3*(n[A]*Pa + n[B]*Pb) + w*P1));
|
||||
|
||||
Faab = k1 * Paab;
|
||||
Fbbc = -k2 * (n[A]*Pabb + n[B]*Pbbb + w*Pbb);
|
||||
Fcca = k3 * (sqr(n[A])*Paaa + 2*n[A]*n[B]*Paab + sqr(n[B])*Pabb
|
||||
+ w*(2*(n[A]*Paa + n[B]*Pab) + w*Pa));
|
||||
}
|
||||
|
||||
|
||||
void Foam::compVolumeIntegrals
|
||||
(
|
||||
const pointField& points,
|
||||
const faceList& faces,
|
||||
const cell& cFaces,
|
||||
const vectorField& fNorm,
|
||||
const scalarField& fW,
|
||||
|
||||
scalar& T0,
|
||||
vector& T1,
|
||||
vector& T2,
|
||||
vector& TP
|
||||
)
|
||||
{
|
||||
T0 = 0;
|
||||
T1 = vector::zero;
|
||||
T2 = vector::zero;
|
||||
TP = vector::zero;
|
||||
|
||||
forAll(cFaces, i)
|
||||
{
|
||||
const vector& n = fNorm[i];
|
||||
|
||||
scalar nx = mag(n[0]);
|
||||
scalar ny = mag(n[1]);
|
||||
scalar nz = mag(n[2]);
|
||||
|
||||
direction A, B, C;
|
||||
|
||||
if (nx > ny && nx > nz)
|
||||
{
|
||||
C = 0;
|
||||
}
|
||||
else
|
||||
{
|
||||
C = (ny > nz) ? 1 : 2;
|
||||
}
|
||||
|
||||
A = (C + 1) % 3;
|
||||
B = (A + 1) % 3;
|
||||
|
||||
scalar Fa, Fb, Fc, Faa, Fbb, Fcc, Faaa, Fbbb, Fccc, Faab, Fbbc, Fcca;
|
||||
compFaceIntegrals
|
||||
(
|
||||
points,
|
||||
faces[cFaces[i]],
|
||||
n,
|
||||
fW[i],
|
||||
A,
|
||||
B,
|
||||
C,
|
||||
|
||||
Fa,
|
||||
Fb,
|
||||
Fc,
|
||||
Faa,
|
||||
Fbb,
|
||||
Fcc,
|
||||
Faaa,
|
||||
Fbbb,
|
||||
Fccc,
|
||||
Faab,
|
||||
Fbbc,
|
||||
Fcca
|
||||
);
|
||||
|
||||
T0 += n[0] * ((A == 0) ? Fa : ((B == 0) ? Fb : Fc));
|
||||
|
||||
T1[A] += n[A] * Faa;
|
||||
T1[B] += n[B] * Fbb;
|
||||
T1[C] += n[C] * Fcc;
|
||||
|
||||
T2[A] += n[A] * Faaa;
|
||||
T2[B] += n[B] * Fbbb;
|
||||
T2[C] += n[C] * Fccc;
|
||||
|
||||
TP[A] += n[A] * Faab;
|
||||
TP[B] += n[B] * Fbbc;
|
||||
TP[C] += n[C] * Fcca;
|
||||
}
|
||||
|
||||
T1 /= 2;
|
||||
T2 /= 3;
|
||||
TP /= 2;
|
||||
}
|
||||
|
||||
|
||||
// Calculate
|
||||
// - r: centre of mass
|
||||
// - J: inertia around origin (point 0,0,0)
|
||||
void Foam::momentOfIntertia
|
||||
(
|
||||
const pointField& points,
|
||||
const faceList& faces,
|
||||
const cell& cFaces,
|
||||
point& r,
|
||||
const pointField& pts,
|
||||
const triFaceList& triFaces,
|
||||
scalar density,
|
||||
scalar& mass,
|
||||
vector& cM,
|
||||
tensor& J
|
||||
)
|
||||
{
|
||||
// Face normals
|
||||
vectorField fNorm(cFaces.size());
|
||||
scalarField fW(cFaces.size());
|
||||
// Reimplemented from: Wm4PolyhedralMassProperties.cpp
|
||||
// File Version: 4.10.0 (2009/11/18)
|
||||
|
||||
forAll(cFaces, i)
|
||||
// Geometric Tools, LC
|
||||
// Copyright (c) 1998-2010
|
||||
// Distributed under the Boost Software License, Version 1.0.
|
||||
// http://www.boost.org/LICENSE_1_0.txt
|
||||
// http://www.geometrictools.com/License/Boost/LICENSE_1_0.txt
|
||||
|
||||
// Boost Software License - Version 1.0 - August 17th, 2003
|
||||
|
||||
// Permission is hereby granted, free of charge, to any person or
|
||||
// organization obtaining a copy of the software and accompanying
|
||||
// documentation covered by this license (the "Software") to use,
|
||||
// reproduce, display, distribute, execute, and transmit the
|
||||
// Software, and to prepare derivative works of the Software, and
|
||||
// to permit third-parties to whom the Software is furnished to do
|
||||
// so, all subject to the following:
|
||||
|
||||
// The copyright notices in the Software and this entire
|
||||
// statement, including the above license grant, this restriction
|
||||
// and the following disclaimer, must be included in all copies of
|
||||
// the Software, in whole or in part, and all derivative works of
|
||||
// the Software, unless such copies or derivative works are solely
|
||||
// in the form of machine-executable object code generated by a
|
||||
// source language processor.
|
||||
|
||||
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
|
||||
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
|
||||
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE AND
|
||||
// NON-INFRINGEMENT. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR
|
||||
// ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE FOR ANY DAMAGES OR
|
||||
// OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
|
||||
// ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE
|
||||
// USE OR OTHER DEALINGS IN THE SOFTWARE.
|
||||
|
||||
const scalar r6 = 1.0/6.0;
|
||||
const scalar r24 = 1.0/24.0;
|
||||
const scalar r60 = 1.0/60.0;
|
||||
const scalar r120 = 1.0/120.0;
|
||||
|
||||
// order: 1, x, y, z, x^2, y^2, z^2, xy, yz, zx
|
||||
scalarField integrals(10, 0.0);
|
||||
|
||||
forAll(triFaces, i)
|
||||
{
|
||||
label faceI = cFaces[i];
|
||||
const triFace& tri(triFaces[i]);
|
||||
|
||||
const face& f = faces[faceI];
|
||||
// vertices of triangle i
|
||||
vector v0 = pts[tri[0]];
|
||||
vector v1 = pts[tri[1]];
|
||||
vector v2 = pts[tri[2]];
|
||||
|
||||
fNorm[i] = f.normal(points);
|
||||
fNorm[i] /= mag(fNorm[i]) + VSMALL;
|
||||
// cross product of edges
|
||||
vector eA = v1 - v0;
|
||||
vector eB = v2 - v0;
|
||||
vector n = eA ^ eB;
|
||||
|
||||
fW[i] = - (fNorm[i] & points[f[0]]);
|
||||
// compute integral terms
|
||||
scalar tmp0, tmp1, tmp2;
|
||||
|
||||
scalar f1x, f2x, f3x, g0x, g1x, g2x;
|
||||
|
||||
tmp0 = v0.x() + v1.x();
|
||||
f1x = tmp0 + v2.x();
|
||||
tmp1 = v0.x()*v0.x();
|
||||
tmp2 = tmp1 + v1.x()*tmp0;
|
||||
f2x = tmp2 + v2.x()*f1x;
|
||||
f3x = v0.x()*tmp1 + v1.x()*tmp2 + v2.x()*f2x;
|
||||
g0x = f2x + v0.x()*(f1x + v0.x());
|
||||
g1x = f2x + v1.x()*(f1x + v1.x());
|
||||
g2x = f2x + v2.x()*(f1x + v2.x());
|
||||
|
||||
scalar f1y, f2y, f3y, g0y, g1y, g2y;
|
||||
|
||||
tmp0 = v0.y() + v1.y();
|
||||
f1y = tmp0 + v2.y();
|
||||
tmp1 = v0.y()*v0.y();
|
||||
tmp2 = tmp1 + v1.y()*tmp0;
|
||||
f2y = tmp2 + v2.y()*f1y;
|
||||
f3y = v0.y()*tmp1 + v1.y()*tmp2 + v2.y()*f2y;
|
||||
g0y = f2y + v0.y()*(f1y + v0.y());
|
||||
g1y = f2y + v1.y()*(f1y + v1.y());
|
||||
g2y = f2y + v2.y()*(f1y + v2.y());
|
||||
|
||||
scalar f1z, f2z, f3z, g0z, g1z, g2z;
|
||||
|
||||
tmp0 = v0.z() + v1.z();
|
||||
f1z = tmp0 + v2.z();
|
||||
tmp1 = v0.z()*v0.z();
|
||||
tmp2 = tmp1 + v1.z()*tmp0;
|
||||
f2z = tmp2 + v2.z()*f1z;
|
||||
f3z = v0.z()*tmp1 + v1.z()*tmp2 + v2.z()*f2z;
|
||||
g0z = f2z + v0.z()*(f1z + v0.z());
|
||||
g1z = f2z + v1.z()*(f1z + v1.z());
|
||||
g2z = f2z + v2.z()*(f1z + v2.z());
|
||||
|
||||
// update integrals
|
||||
integrals[0] += n.x()*f1x;
|
||||
integrals[1] += n.x()*f2x;
|
||||
integrals[2] += n.y()*f2y;
|
||||
integrals[3] += n.z()*f2z;
|
||||
integrals[4] += n.x()*f3x;
|
||||
integrals[5] += n.y()*f3y;
|
||||
integrals[6] += n.z()*f3z;
|
||||
integrals[7] += n.x()*(v0.y()*g0x + v1.y()*g1x + v2.y()*g2x);
|
||||
integrals[8] += n.y()*(v0.z()*g0y + v1.z()*g1y + v2.z()*g2y);
|
||||
integrals[9] += n.z()*(v0.x()*g0z + v1.x()*g1z + v2.x()*g2z);
|
||||
}
|
||||
|
||||
integrals[0] *= r6;
|
||||
integrals[1] *= r24;
|
||||
integrals[2] *= r24;
|
||||
integrals[3] *= r24;
|
||||
integrals[4] *= r60;
|
||||
integrals[5] *= r60;
|
||||
integrals[6] *= r60;
|
||||
integrals[7] *= r120;
|
||||
integrals[8] *= r120;
|
||||
integrals[9] *= r120;
|
||||
|
||||
scalar T0;
|
||||
vector T1, T2, TP;
|
||||
// mass
|
||||
mass = integrals[0];
|
||||
|
||||
compVolumeIntegrals
|
||||
(
|
||||
points,
|
||||
faces,
|
||||
cFaces,
|
||||
fNorm,
|
||||
fW,
|
||||
// center of mass
|
||||
cM = vector(integrals[1], integrals[2], integrals[3])/mass;
|
||||
|
||||
T0,
|
||||
T1,
|
||||
T2,
|
||||
TP
|
||||
);
|
||||
// inertia relative to origin
|
||||
J.xx() = integrals[5] + integrals[6];
|
||||
J.xy() = -integrals[7];
|
||||
J.xz() = -integrals[9];
|
||||
J.yx() = J.xy();
|
||||
J.yy() = integrals[4] + integrals[6];
|
||||
J.yz() = -integrals[8];
|
||||
J.zx() = J.xz();
|
||||
J.zy() = J.yz();
|
||||
J.zz() = integrals[4] + integrals[5];
|
||||
|
||||
const scalar density = 1.0; /* assume unit density */
|
||||
// inertia relative to center of mass
|
||||
J -= mass*((cM & cM)*I - cM*cM);
|
||||
|
||||
scalar mass = density * T0;
|
||||
|
||||
/* compute center of mass */
|
||||
r = T1 / T0;
|
||||
|
||||
/* compute inertia tensor */
|
||||
J.xx() = density * (T2[1] + T2[2]);
|
||||
J.yy() = density * (T2[2] + T2[0]);
|
||||
J.zz() = density * (T2[0] + T2[1]);
|
||||
J.xy() = J.yx() = - density * TP[0];
|
||||
J.yz() = J.zy() = - density * TP[1];
|
||||
J.zx() = J.xz() = - density * TP[2];
|
||||
|
||||
///* translate inertia tensor to center of mass */
|
||||
//J[XX] -= mass * (r[1]*r[1] + r[2]*r[2]);
|
||||
//J[YY] -= mass * (r[2]*r[2] + r[0]*r[0]);
|
||||
//J[ZZ] -= mass * (r[0]*r[0] + r[1]*r[1]);
|
||||
//J[XY] = J[YX] += mass * r[0] * r[1];
|
||||
//J[YZ] = J[ZY] += mass * r[1] * r[2];
|
||||
//J[ZX] = J[XZ] += mass * r[2] * r[0];
|
||||
// Apply density
|
||||
mass *= density;
|
||||
J *= density;
|
||||
}
|
||||
|
||||
|
||||
void Foam::momentOfInertia::massPropertiesShell
|
||||
(
|
||||
const pointField& pts,
|
||||
const triFaceList& triFaces,
|
||||
scalar density,
|
||||
scalar& mass,
|
||||
vector& cM,
|
||||
tensor& J
|
||||
)
|
||||
{
|
||||
// Reset properties for accumulation
|
||||
|
||||
mass = 0.0;
|
||||
cM = vector::zero;
|
||||
J = tensor::zero;
|
||||
|
||||
// Find centre of mass
|
||||
|
||||
forAll(triFaces, i)
|
||||
{
|
||||
const triFace& tri(triFaces[i]);
|
||||
|
||||
triPointRef t
|
||||
(
|
||||
pts[tri[0]],
|
||||
pts[tri[1]],
|
||||
pts[tri[2]]
|
||||
);
|
||||
|
||||
scalar triMag = t.mag();
|
||||
|
||||
cM += triMag*t.centre();
|
||||
|
||||
mass += triMag;
|
||||
}
|
||||
|
||||
cM /= mass;
|
||||
|
||||
mass *= density;
|
||||
|
||||
// Find inertia around centre of mass
|
||||
|
||||
forAll(triFaces, i)
|
||||
{
|
||||
const triFace& tri(triFaces[i]);
|
||||
|
||||
J += triPointRef
|
||||
(
|
||||
pts[tri[0]],
|
||||
pts[tri[1]],
|
||||
pts[tri[2]]
|
||||
).inertia(cM, density);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void Foam::momentOfInertia::massPropertiesSolid
|
||||
(
|
||||
const triSurface& surf,
|
||||
scalar density,
|
||||
scalar& mass,
|
||||
vector& cM,
|
||||
tensor& J
|
||||
)
|
||||
{
|
||||
triFaceList faces(surf.size());
|
||||
|
||||
forAll(surf, i)
|
||||
{
|
||||
faces[i] = triFace(surf[i]);
|
||||
}
|
||||
|
||||
massPropertiesSolid(surf.points(), faces, density, mass, cM, J);
|
||||
}
|
||||
|
||||
|
||||
void Foam::momentOfInertia::massPropertiesShell
|
||||
(
|
||||
const triSurface& surf,
|
||||
scalar density,
|
||||
scalar& mass,
|
||||
vector& cM,
|
||||
tensor& J
|
||||
)
|
||||
{
|
||||
triFaceList faces(surf.size());
|
||||
|
||||
forAll(surf, i)
|
||||
{
|
||||
faces[i] = triFace(surf[i]);
|
||||
}
|
||||
|
||||
massPropertiesShell(surf.points(), faces, density, mass, cM, J);
|
||||
}
|
||||
|
||||
|
||||
Foam::tensor Foam::momentOfInertia::applyParallelAxisTheorem
|
||||
(
|
||||
scalar mass,
|
||||
const vector& cM,
|
||||
const tensor& J,
|
||||
const vector& refPt
|
||||
)
|
||||
{
|
||||
// The displacement vector (refPt = cM) is the displacement of the
|
||||
// new reference point from the centre of mass of the body that
|
||||
// the inertia tensor applies to.
|
||||
|
||||
vector d = (refPt - cM);
|
||||
|
||||
return J + mass*((d & d)*I - d*d);
|
||||
}
|
||||
|
||||
|
||||
Foam::tmp<Foam::tensorField> Foam::momentOfInertia::meshInertia
|
||||
(
|
||||
const polyMesh& mesh
|
||||
)
|
||||
{
|
||||
tmp<tensorField> tTf = tmp<tensorField>(new tensorField(mesh.nCells()));
|
||||
|
||||
tensorField& tf = tTf();
|
||||
|
||||
forAll(tf, cI)
|
||||
{
|
||||
tf[cI] = meshInertia(mesh, cI);
|
||||
}
|
||||
|
||||
return tTf;
|
||||
}
|
||||
|
||||
|
||||
Foam::tensor Foam::momentOfInertia::meshInertia
|
||||
(
|
||||
const polyMesh& mesh,
|
||||
label cellI
|
||||
)
|
||||
{
|
||||
List<tetIndices> cellTets = polyMeshTetDecomposition::cellTetIndices
|
||||
(
|
||||
mesh,
|
||||
cellI
|
||||
);
|
||||
|
||||
triFaceList faces(cellTets.size());
|
||||
|
||||
forAll(cellTets, cTI)
|
||||
{
|
||||
faces[cTI] = cellTets[cTI].faceTriIs(mesh);
|
||||
}
|
||||
|
||||
scalar m = 0.0;
|
||||
vector cM = vector::zero;
|
||||
tensor J = tensor::zero;
|
||||
|
||||
massPropertiesSolid(mesh.points(), faces, 1.0, m, cM, J);
|
||||
|
||||
return J;
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
|
||||
@ -25,6 +25,9 @@ Class
|
||||
momentOfInertia
|
||||
|
||||
Description
|
||||
Calculates the inertia tensor and principal axes and moments of a
|
||||
polyhedra/cells/triSurfaces. Inertia can either be of the solid body or
|
||||
of a thin shell.
|
||||
|
||||
SourceFiles
|
||||
momentOfInertia.H
|
||||
@ -34,34 +37,86 @@ SourceFiles
|
||||
#ifndef momentOfInertia_H
|
||||
#define momentOfInertia_H
|
||||
|
||||
#include "tensor.H"
|
||||
#include "primitiveMesh.H"
|
||||
#include "tetPointRef.H"
|
||||
#include "triFaceList.H"
|
||||
#include "triSurface.H"
|
||||
#include "polyMesh.H"
|
||||
#include "polyMeshTetDecomposition.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
namespace Foam
|
||||
{
|
||||
|
||||
////- Moment of inertia around cell centre for single cell.
|
||||
//tensor momentOfInertia
|
||||
//(
|
||||
// const pointField&,
|
||||
// const faceList&,
|
||||
// const cell&,
|
||||
// const point& cc
|
||||
//);
|
||||
/*---------------------------------------------------------------------------*\
|
||||
Class momentOfInertia Declaration
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
class momentOfInertia
|
||||
{
|
||||
|
||||
public:
|
||||
|
||||
static void massPropertiesSolid
|
||||
(
|
||||
const pointField& pts,
|
||||
const triFaceList& triFaces,
|
||||
scalar density,
|
||||
scalar& mass,
|
||||
vector& cM,
|
||||
tensor& J
|
||||
);
|
||||
|
||||
static void massPropertiesShell
|
||||
(
|
||||
const pointField& pts,
|
||||
const triFaceList& triFaces,
|
||||
scalar density,
|
||||
scalar& mass,
|
||||
vector& cM,
|
||||
tensor& J
|
||||
);
|
||||
|
||||
static void massPropertiesSolid
|
||||
(
|
||||
const triSurface& surf,
|
||||
scalar density,
|
||||
scalar& mass,
|
||||
vector& cM,
|
||||
tensor& J
|
||||
);
|
||||
|
||||
static void massPropertiesShell
|
||||
(
|
||||
const triSurface& surf,
|
||||
scalar density,
|
||||
scalar& mass,
|
||||
vector& cM,
|
||||
tensor& J
|
||||
);
|
||||
|
||||
static tensor applyParallelAxisTheorem
|
||||
(
|
||||
scalar mass,
|
||||
const vector& cM,
|
||||
const tensor& J,
|
||||
const vector& refPt
|
||||
);
|
||||
|
||||
// Calculate the inertia tensor for all cells in the mesh
|
||||
static tmp<tensorField> meshInertia
|
||||
(
|
||||
const polyMesh& mesh
|
||||
);
|
||||
|
||||
// Calculate the inertia tensor the given cell
|
||||
static tensor meshInertia
|
||||
(
|
||||
const polyMesh& mesh,
|
||||
label cellI
|
||||
);
|
||||
};
|
||||
|
||||
// Calculate
|
||||
// - centre of mass
|
||||
// - inertia tensor around (0,0,0)
|
||||
void momentOfIntertia
|
||||
(
|
||||
const pointField&,
|
||||
const faceList&,
|
||||
const cell&,
|
||||
point& r,
|
||||
tensor& Jorigin
|
||||
);
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
|
||||
@ -34,7 +34,7 @@ divSchemes
|
||||
div(phi,epsilon) Gauss upwind;
|
||||
div(phi,R) Gauss upwind;
|
||||
div(R) Gauss linear;
|
||||
div((nuEff*dev(grad(U).T()))) Gauss linear;
|
||||
div((nuEff*dev(T(grad(U))))) Gauss linear;
|
||||
}
|
||||
|
||||
laplacianSchemes
|
||||
|
||||
@ -32,7 +32,7 @@ divSchemes
|
||||
div(phi,T) Gauss upwind;
|
||||
div(phi,k) Gauss upwind;
|
||||
div(phi,epsilon) Gauss upwind;
|
||||
div((nuEff*dev(grad(U).T()))) Gauss linear;
|
||||
div((nuEff*dev(T(grad(U))))) Gauss linear;
|
||||
}
|
||||
|
||||
laplacianSchemes
|
||||
|
||||
@ -34,7 +34,7 @@ divSchemes
|
||||
div(phi,epsilon) Gauss upwind;
|
||||
div(phi,R) Gauss upwind;
|
||||
div(R) Gauss linear;
|
||||
div((nuEff*dev(grad(U).T()))) Gauss linear;
|
||||
div((nuEff*dev(T(grad(U))))) Gauss linear;
|
||||
}
|
||||
|
||||
laplacianSchemes
|
||||
|
||||
@ -33,7 +33,7 @@ divSchemes
|
||||
div(phi,U) Gauss limitedLinearV 1;
|
||||
div(phi,k) Gauss limitedLinear 1;
|
||||
div(phi,epsilon) Gauss limitedLinear 1;
|
||||
div((nuEff*dev(grad(U).T()))) Gauss linear;
|
||||
div((nuEff*dev(T(grad(U))))) Gauss linear;
|
||||
}
|
||||
|
||||
laplacianSchemes
|
||||
|
||||
@ -35,7 +35,7 @@ divSchemes
|
||||
div(phi,U) Gauss upwind;
|
||||
div(phi,k) Gauss upwind;
|
||||
div(phi,epsilon) Gauss upwind;
|
||||
div((nuEff*dev(grad(U).T()))) Gauss linear;
|
||||
div((nuEff*dev(T(grad(U))))) Gauss linear;
|
||||
|
||||
div(-phi,Ua) Gauss upwind;
|
||||
div((nuEff*dev(grad(Ua).T()))) Gauss linear;
|
||||
|
||||
@ -32,7 +32,7 @@ divSchemes
|
||||
div(phi,epsilon) Gauss linear;
|
||||
div(phi,R) Gauss linear;
|
||||
div(phi,nuTilda) Gauss linear;
|
||||
div((nuEff*dev(grad(U).T()))) Gauss linear;
|
||||
div((nuEff*dev(T(grad(U))))) Gauss linear;
|
||||
}
|
||||
|
||||
laplacianSchemes
|
||||
|
||||
@ -32,7 +32,7 @@ divSchemes
|
||||
div(phi,epsilon) Gauss linear;
|
||||
div(phi,R) Gauss linear;
|
||||
div(phi,nuTilda) Gauss linear;
|
||||
div((nuEff*dev(grad(U).T()))) Gauss linear;
|
||||
div((nuEff*dev(T(grad(U))))) Gauss linear;
|
||||
}
|
||||
|
||||
laplacianSchemes
|
||||
|
||||
@ -33,7 +33,7 @@ divSchemes
|
||||
div(phi,omega) Gauss linear;
|
||||
div(phi,R) Gauss linear;
|
||||
div(phi,nuTilda) Gauss linear;
|
||||
div((nuEff*dev(grad(U).T()))) Gauss linear;
|
||||
div((nuEff*dev(T(grad(U))))) Gauss linear;
|
||||
}
|
||||
|
||||
laplacianSchemes
|
||||
|
||||
@ -35,7 +35,7 @@ divSchemes
|
||||
div(phi,B) Gauss limitedLinear 1;
|
||||
div(B) Gauss linear;
|
||||
div(phi,nuTilda) Gauss limitedLinear 1;
|
||||
div((nuEff*dev(grad(U).T()))) Gauss linear;
|
||||
div((nuEff*dev(T(grad(U))))) Gauss linear;
|
||||
}
|
||||
|
||||
laplacianSchemes
|
||||
|
||||
@ -30,7 +30,7 @@ divSchemes
|
||||
{
|
||||
default none;
|
||||
div(phi,U) Gauss linear;
|
||||
div((nuEff*dev(grad(U).T()))) Gauss linear;
|
||||
div((nuEff*dev(T(grad(U))))) Gauss linear;
|
||||
}
|
||||
|
||||
laplacianSchemes
|
||||
|
||||
@ -32,7 +32,7 @@ divSchemes
|
||||
div(phi,U) Gauss linearUpwind grad(U);
|
||||
div(phi,k) Gauss limitedLinear 1;
|
||||
div(phi,omega) Gauss limitedLinear 1;
|
||||
div((nuEff*dev(grad(U).T()))) Gauss linear;
|
||||
div((nuEff*dev(T(grad(U))))) Gauss linear;
|
||||
}
|
||||
|
||||
laplacianSchemes
|
||||
|
||||
@ -32,7 +32,7 @@ divSchemes
|
||||
div(phi,U) Gauss linearUpwind grad(U);
|
||||
div(phi,k) Gauss upwind;
|
||||
div(phi,omega) Gauss upwind;
|
||||
div((nuEff*dev(grad(U).T()))) Gauss linear;
|
||||
div((nuEff*dev(T(grad(U))))) Gauss linear;
|
||||
}
|
||||
|
||||
laplacianSchemes
|
||||
|
||||
@ -36,7 +36,7 @@ divSchemes
|
||||
div(phi,R) Gauss limitedLinear 1;
|
||||
div(R) Gauss linear;
|
||||
div(phi,nuTilda) Gauss limitedLinear 1;
|
||||
div((nuEff*dev(grad(U).T()))) Gauss linear;
|
||||
div((nuEff*dev(T(grad(U))))) Gauss linear;
|
||||
}
|
||||
|
||||
laplacianSchemes
|
||||
|
||||
@ -36,7 +36,7 @@ divSchemes
|
||||
div(phi,R) Gauss limitedLinear 1;
|
||||
div(R) Gauss linear;
|
||||
div(phi,nuTilda) Gauss limitedLinear 1;
|
||||
div((nuEff*dev(grad(U).T()))) Gauss linear;
|
||||
div((nuEff*dev(T(grad(U))))) Gauss linear;
|
||||
}
|
||||
|
||||
laplacianSchemes
|
||||
|
||||
@ -35,7 +35,7 @@ divSchemes
|
||||
div(phi,B) Gauss limitedLinear 1;
|
||||
div(phi,nuTilda) Gauss limitedLinear 1;
|
||||
div(B) Gauss linear;
|
||||
div((nuEff*dev(grad(U).T()))) Gauss linear;
|
||||
div((nuEff*dev(T(grad(U))))) Gauss linear;
|
||||
}
|
||||
|
||||
laplacianSchemes
|
||||
|
||||
@ -35,7 +35,7 @@ divSchemes
|
||||
div(phi,B) Gauss limitedLinear 1;
|
||||
div(phi,nuTilda) Gauss limitedLinear 1;
|
||||
div(B) Gauss linear;
|
||||
div((nuEff*dev(grad(U).T()))) Gauss linear;
|
||||
div((nuEff*dev(T(grad(U))))) Gauss linear;
|
||||
}
|
||||
|
||||
laplacianSchemes
|
||||
|
||||
@ -36,7 +36,7 @@ divSchemes
|
||||
div(phi,R) Gauss limitedLinear 1;
|
||||
div(R) Gauss linear;
|
||||
div(phi,nuTilda) Gauss limitedLinear 1;
|
||||
div((nuEff*dev(grad(U).T()))) Gauss linear;
|
||||
div((nuEff*dev(T(grad(U))))) Gauss linear;
|
||||
}
|
||||
|
||||
laplacianSchemes
|
||||
|
||||
@ -32,7 +32,7 @@ divSchemes
|
||||
default none;
|
||||
div(phi,U) Gauss linearUpwind Gauss linear;
|
||||
div(phi,nuTilda) Gauss linearUpwind Gauss linear;
|
||||
div((nuEff*dev(grad(U).T()))) Gauss linear;
|
||||
div((nuEff*dev(T(grad(U))))) Gauss linear;
|
||||
}
|
||||
|
||||
laplacianSchemes
|
||||
|
||||
@ -32,7 +32,7 @@ divSchemes
|
||||
div(phi,U) Gauss linearUpwind grad(U);
|
||||
div(phi,k) Gauss upwind;
|
||||
div(phi,omega) Gauss upwind;
|
||||
div((nuEff*dev(grad(U).T()))) Gauss linear;
|
||||
div((nuEff*dev(T(grad(U))))) Gauss linear;
|
||||
}
|
||||
|
||||
laplacianSchemes
|
||||
|
||||
@ -19,4 +19,7 @@ rm -rf constant/polyMesh/sets
|
||||
#rm -rf constant/refinementHistory
|
||||
#rm -rf constant/surfaceIndex
|
||||
|
||||
# Reset decomposeParDict
|
||||
cp system/decomposeParDict-nonPar system/decomposeParDict
|
||||
|
||||
# ----------------------------------------------------------------- end-of-file
|
||||
|
||||
@ -7,6 +7,7 @@ cd ${0%/*} || exit 1 # run from this directory
|
||||
compileApplication windSimpleFoam
|
||||
|
||||
runApplication blockMesh
|
||||
cp system/decomposeParDict-nonPar system/decomposeParDict
|
||||
runApplication decomposePar
|
||||
|
||||
#runApplication snappyHexMesh -overwrite
|
||||
@ -14,6 +15,7 @@ runApplication decomposePar
|
||||
#runApplication setsToZones -noFlipMap
|
||||
#runApplication windSimpleFoam
|
||||
|
||||
cp system/decomposeParDict-par system/decomposeParDict
|
||||
runParallel snappyHexMesh 2 -overwrite
|
||||
|
||||
# Add wildcard entries for meshes patches since not preserved
|
||||
|
||||
@ -0,0 +1,29 @@
|
||||
/*--------------------------------*- C++ -*----------------------------------*\
|
||||
| ========= | |
|
||||
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
|
||||
| \\ / O peration | Version: dev |
|
||||
| \\ / A nd | Web: www.OpenFOAM.com |
|
||||
| \\/ M anipulation | |
|
||||
\*---------------------------------------------------------------------------*/
|
||||
FoamFile
|
||||
{
|
||||
version 2.0;
|
||||
format ascii;
|
||||
class dictionary;
|
||||
object decomposeParDict;
|
||||
}
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
numberOfSubdomains 2;
|
||||
|
||||
method hierarchical;
|
||||
|
||||
hierarchicalCoeffs
|
||||
{
|
||||
n (2 1 1);
|
||||
delta 0.001;
|
||||
order xyz;
|
||||
}
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -0,0 +1,22 @@
|
||||
/*--------------------------------*- C++ -*----------------------------------*\
|
||||
| ========= | |
|
||||
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
|
||||
| \\ / O peration | Version: dev |
|
||||
| \\ / A nd | Web: www.OpenFOAM.com |
|
||||
| \\/ M anipulation | |
|
||||
\*---------------------------------------------------------------------------*/
|
||||
FoamFile
|
||||
{
|
||||
version 2.0;
|
||||
format ascii;
|
||||
class dictionary;
|
||||
object decomposeParDict;
|
||||
}
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
numberOfSubdomains 2;
|
||||
|
||||
method ptscotch;
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -30,7 +30,7 @@ divSchemes
|
||||
{
|
||||
default none;
|
||||
div(phi,U) Gauss upwind grad(U);
|
||||
div((nuEff*dev(grad(U).T()))) Gauss linear;
|
||||
div((nuEff*dev(T(grad(U))))) Gauss linear;
|
||||
div(phi,epsilon) Gauss upwind;
|
||||
div(phi,k) Gauss upwind;
|
||||
}
|
||||
|
||||
@ -32,7 +32,7 @@ divSchemes
|
||||
div(phirb,alpha) Gauss interfaceCompression 1;
|
||||
div(phi,p_rgh) Gauss upwind;
|
||||
div(phi,k) Gauss vanLeer;
|
||||
div((nuEff*dev(grad(U).T()))) Gauss linear;
|
||||
div((nuEff*dev(T(grad(U))))) Gauss linear;
|
||||
}
|
||||
|
||||
laplacianSchemes
|
||||
|
||||
@ -32,7 +32,7 @@ divSchemes
|
||||
div(phirb,alpha) Gauss interfaceCompression 1;
|
||||
div(phi,p_rgh) Gauss upwind;
|
||||
div(phi,k) Gauss vanLeer;
|
||||
div((nuEff*dev(grad(U).T()))) Gauss linear;
|
||||
div((nuEff*dev(T(grad(U))))) Gauss linear;
|
||||
}
|
||||
|
||||
laplacianSchemes
|
||||
|
||||
@ -34,7 +34,7 @@ divSchemes
|
||||
div(phi,B) Gauss limitedLinear 1;
|
||||
div(B) Gauss linear;
|
||||
div(phi,nuTilda) Gauss limitedLinear 1;
|
||||
div((nuEff*dev(grad(U).T()))) Gauss linear;
|
||||
div((nuEff*dev(T(grad(U))))) Gauss linear;
|
||||
}
|
||||
|
||||
laplacianSchemes
|
||||
|
||||
@ -35,7 +35,7 @@ divSchemes
|
||||
div(phi,R) Gauss upwind;
|
||||
div(R) Gauss linear;
|
||||
div(phi,nuTilda) Gauss upwind;
|
||||
div((nuEff*dev(grad(U).T()))) Gauss linear;
|
||||
div((nuEff*dev(T(grad(U))))) Gauss linear;
|
||||
}
|
||||
|
||||
laplacianSchemes
|
||||
|
||||
Reference in New Issue
Block a user