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;
101 pointNormals[pI].normalise();
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);
void swap(UList< T > &list) noexcept
Swap content with another UList of the same type in constant time.
label nPoints() const
Number of points supporting patch faces.
void size(const label n)
Older name for setAddressableSize.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
T & ref() const
Return non-const reference to the contents of a non-null managed pointer.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Vector2D< Cmpt > x() const
Extract vector for row 0.
Vector2D< scalar > vector2D
A 2D vector of scalars obtained from the generic Vector2D.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Ignore writing from objectRegistry::writeObject()
Geometric class that creates a 3D plane and can return the intersection point between a line and the ...
quaternion normalised(const quaternion &q)
Return the normalised (unit) quaternion of the given quaternion.
tensor rotationTensor(const vector &n1, const vector &n2)
Rotational transformation tensor from vector n1 to n2.
Vector2D< Cmpt > y() const
Extract vector for row 1.
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
A templated (2 x 2) tensor of objects of <T> derived from VectorSpace.
const labelList & meshPoints() const
Return labelList of mesh points in patch.
bool set(const direction d) const
Is the vector in the direction d set.
#define forAll(list, i)
Loop across all elements in list.
void normalize()
Same as normalise.
const Map< label > & meshPointMap() const
Mesh point map.
const dimensionedScalar e
Elementary charge.
An edge is a list of two vertex labels. This can correspond to a directed graph edge or an edge on a ...
A class for handling words, derived from Foam::string.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
A triangular face using a FixedList of labels corresponding to mesh vertices.
static tmp< T > New(Args &&... args)
Construct tmp with forwarding arguments.
const Field< point_type > & points() const noexcept
Return reference to global points.
errorManip< error > abort(error &err)
A 1D vector of objects of type <T>, where the size of the vector is known and can be used for subscri...
const labelListList & pointFaces() const
Return point-face addressing.
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
const word & constant() const noexcept
Return constant name.
void T(FieldField< Field, Type > &f1, const FieldField< Field, Type > &f2)
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]
A templated (N x N) square matrix of objects of <Type>, containing N*N elements, derived from Matrix...
Representation of a 3D Cartesian coordinate system as a Vector of row vectors.
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
A templated (2 x 2) symmetric tensor of objects of <T>, effectively containing 3 elements, derived from VectorSpace.
messageStream Info
Information stream (stdout output on master, null elsewhere)
dimensionedVector eigenValues(const dimensionedSymmTensor &dt)
void clear() const noexcept
If object pointer points to valid object: delete object and set pointer to nullptr.
A class for managing temporary objects.
Triangulated surface description with patch information.
Tensor of scalars, i.e. Tensor<scalar>.
Defines the attributes of an object for which implicit objectRegistry management is supported...
static const SymmTensor2D< Cmpt > zero
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
void LUsolve(scalarSquareMatrix &matrix, List< Type > &source)
Solve the matrix using LU decomposition with pivoting returning the LU form of the matrix and the sol...
static constexpr const zero Zero
Global zero (0)