50 label index =
f.
find(pI);
71 Info<<
"Calculating vertex normals" <<
endl;
74 auto& pointNormals = tpointNormals.ref();
84 for (
const label facei :
pFaces)
90 const scalar weight = vertexNormalWeight
98 pointNormals[pI] += weight * areaNorm;
104 return tpointNormals;
119 auto& pointTriads = tpointTriads.ref();
124 const vector& normal = pointNormals[meshPointMap[pI]];
126 if (
mag(normal) < SMALL)
138 pointTriads[meshPointMap[pI]] =
triad(dir1, dir2, normal);
153 Info<<
"Calculating face curvature" <<
endl;
178 const edge&
e = fEdges[feI];
180 edgeVectors[feI] =
e.vec(
points);
181 normalDifferences[feI] =
182 pointNormals[meshPointMap[
e[0]]]
183 - pointNormals[meshPointMap[
e[1]]];
187 const vector& e0 = edgeVectors[0];
189 const vector e1 = (e0 ^ eN);
191 if (
magSqr(eN) < ROOTVSMALL)
196 triad faceCoordSys(e0, e1, eN);
204 for (label i = 0; i < 3; ++i)
206 scalar
x = edgeVectors[i] & faceCoordSys[0];
207 scalar
y = edgeVectors[i] & faceCoordSys[1];
215 scalar dndx = normalDifferences[i] & faceCoordSys[0];
216 scalar dndy = normalDifferences[i] & faceCoordSys[1];
219 Z[1] += dndx*
y + dndy*
x;
232 const label patchPointIndex = meshPointMap[
f[fpI]];
234 const triad& ptCoordSys = pointTriads[patchPointIndex];
236 if (!ptCoordSys.
set())
243 triad rotatedFaceCoordSys = rotTensor &
tensor(faceCoordSys);
249 ptCoordSys[0] & rotatedFaceCoordSys[0],
250 ptCoordSys[0] & rotatedFaceCoordSys[1]
254 ptCoordSys[1] & rotatedFaceCoordSys[0],
255 ptCoordSys[1] & rotatedFaceCoordSys[1]
266 projTensor.
x() & (secondFundamentalTensor & projTensor.
x()),
267 projTensor.
x() & (secondFundamentalTensor & projTensor.
y()),
268 projTensor.
y() & (secondFundamentalTensor & projTensor.
y())
276 meshPoints[patchPointIndex],
282 pointFundamentalTensors[patchPointIndex] +=
283 weight*projectedFundamentalTensor;
285 accumulatedWeights[patchPointIndex] += weight;
293 Info<<
"edgeVecs = " << edgeVectors[0] <<
" "
294 << edgeVectors[1] <<
" "
295 << edgeVectors[2] <<
endl;
296 Info<<
"normDiff = " << normalDifferences[0] <<
" "
297 << normalDifferences[1] <<
" "
298 << normalDifferences[2] <<
endl;
299 Info<<
"faceCoordSys = " << faceCoordSys <<
endl;
306 scalarField& curvatureAtPoints = tcurvatureAtPoints.ref();
308 forAll(curvatureAtPoints, pI)
310 pointFundamentalTensors[pI] /= (accumulatedWeights[pI] + SMALL);
316 scalar curvature =
max
318 mag(principalCurvatures[0]),
319 mag(principalCurvatures[1])
323 curvatureAtPoints[meshPoints[pI]] = curvature;
326 return tcurvatureAtPoints;
351 const word& basename,
355 Info<<
"Extracting curvature of surface at the points." <<
endl;
364 basename +
".curvature",
376 outputField.
swap(curv);
378 outputField.
swap(curv);
Defines the attributes of an object for which implicit objectRegistry management is supported,...
A HashTable to objects of type <T> with a label key.
label nPoints() const
Number of points supporting patch faces.
const Map< label > & meshPointMap() const
Mesh point map.
const labelList & meshPoints() const
Return labelList of mesh points in patch.
const Field< point_type > & points() const noexcept
Return reference to global points.
const labelListList & pointFaces() const
Return point-face addressing.
A templated (2 x 2) symmetric tensor of objects of <T>, effectively containing 3 elements,...
A templated (N x N) square matrix of objects of <Type>, containing N*N elements, derived from Matrix.
A templated (2 x 2) tensor of objects of <T> derived from VectorSpace.
Vector2D< Cmpt > y() const
Extract vector for row 1.
Vector2D< Cmpt > x() const
Extract vector for row 0.
const word & constant() const
Return constant name.
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
static autoPtr< Time > New()
Construct (dummy) Time - no functionObjects or libraries.
A 1D vector of objects of type <T>, where the size of the vector is known and can be used for subscri...
void swap(UList< T > &list)
Swap content with another UList of the same type in constant time.
label find(const T &val, label pos=0) const
Find index of the first occurrence of the value.
label rcIndex(const label i) const noexcept
void size(const label n)
Older name for setAddressableSize.
label fcIndex(const label i) const noexcept
Vector< Cmpt > & normalise(const scalar tol=ROOTVSMALL)
Inplace normalise the vector by its magnitude.
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Geometric class that creates a 3D plane and can return the intersection point between a line and the ...
virtual bool write(const bool valid=true) const
Write using setting from DB.
A class for managing temporary objects.
void clear() const noexcept
A triangular face using a FixedList of labels corresponding to mesh vertices.
Triangulated surface description with patch information.
Representation of a 3D Cartesian coordinate system as a Vector of row vectors.
bool set(const direction d) const
Is the vector in the direction d set.
void normalize()
Same as normalise.
A class for handling words, derived from Foam::string.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
quaternion normalised(const quaternion &q)
Return the normalised (unit) quaternion of the given quaternion.
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
dimensionedVector eigenValues(const dimensionedSymmTensor &dt)
tensor rotationTensor(const vector &n1, const vector &n2)
Rotational transformation tensor from vector n1 to n2.
messageStream Info
Information stream (stdout output on master, null elsewhere)
Ostream & endl(Ostream &os)
Add newline and flush stream.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
errorManip< error > abort(error &err)
static constexpr const zero Zero
Global zero (0)
void LUsolve(scalarSquareMatrix &matrix, List< Type > &source)
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
Info<< "Finished reading KIVA file"<< endl;cellShapeList cellShapes(nPoints);labelList cellZoning(nPoints, -1);const cellModel &hex=cellModel::ref(cellModel::HEX);labelList hexLabels(8);label activeCells=0;labelList pointMap(nPoints);forAll(pointMap, i){ pointMap[i]=i;}for(label i=0;i< nPoints;i++){ if(f[i] > 0.0) { hexLabels[0]=i;hexLabels[1]=i1tab[i];hexLabels[2]=i3tab[i1tab[i]];hexLabels[3]=i3tab[i];hexLabels[4]=i8tab[i];hexLabels[5]=i1tab[i8tab[i]];hexLabels[6]=i3tab[i1tab[i8tab[i]]];hexLabels[7]=i3tab[i8tab[i]];cellShapes[activeCells].reset(hex, hexLabels);edgeList edges=cellShapes[activeCells].edges();forAll(edges, ei) { if(edges[ei].mag(points)< SMALL) { label start=pointMap[edges[ei].start()];while(start !=pointMap[start]) { start=pointMap[start];} label end=pointMap[edges[ei].end()];while(end !=pointMap[end]) { end=pointMap[end];} label minLabel=min(start, end);pointMap[start]=pointMap[end]=minLabel;} } cellZoning[activeCells]=idreg[i];activeCells++;}}cellShapes.setSize(activeCells);cellZoning.setSize(activeCells);forAll(cellShapes, celli){ cellShape &cs=cellShapes[celli];forAll(cs, i) { cs[i]=pointMap[cs[i]];} cs.collapse();}label bcIDs[11]={-1, 0, 2, 4, -1, 5, -1, 6, 7, 8, 9};const label nBCs=12;const word *kivaPatchTypes[nBCs]={ &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &symmetryPolyPatch::typeName, &wedgePolyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &symmetryPolyPatch::typeName, &oldCyclicPolyPatch::typeName};enum patchTypeNames{ PISTON, VALVE, LINER, CYLINDERHEAD, AXIS, WEDGE, INFLOW, OUTFLOW, PRESIN, PRESOUT, SYMMETRYPLANE, CYCLIC};const char *kivaPatchNames[nBCs]={ "piston", "valve", "liner", "cylinderHead", "axis", "wedge", "inflow", "outflow", "presin", "presout", "symmetryPlane", "cyclic"};List< SLList< face > > pFaces[nBCs]
#define forAll(list, i)
Loop across all elements in list.
A non-counting (dummy) refCount.