Go to the documentation of this file.
50 label index =
f.find(pI);
71 Info<<
"Calculating vertex normals" <<
endl;
74 auto& pointNormals = tpointNormals.ref();
84 for (
const label facei :
pFaces)
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);
const Field< point_type > & points() const noexcept
Return reference to global points.
void LUsolve(scalarSquareMatrix &matrix, List< Type > &source)
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Defines the attributes of an object for which implicit objectRegistry management is supported,...
const labelListList & pointFaces() const
Return point-face addressing.
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
A class for handling words, derived from Foam::string.
void clear() const noexcept
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
A class for managing temporary objects.
static constexpr const zero Zero
Global zero (0)
dimensionedVector eigenValues(const dimensionedSymmTensor &dt)
A templated (2 x 2) symmetric tensor of objects of <T>, effectively containing 3 elements,...
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
void normalize()
Normalize each set axis vector to have a unit magnitude.
Ostream & endl(Ostream &os)
Add newline and flush stream.
#define forAll(list, i)
Loop across all elements in list.
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
A templated (2 x 2) tensor of objects of <T> derived from VectorSpace.
Geometric class that creates a 3D plane and can return the intersection point between a line and the ...
Triangulated surface description with patch information.
messageStream Info
Information stream (stdout output on master, null elsewhere)
void T(FieldField< Field, Type > &f1, const FieldField< Field, Type > &f2)
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
label nPoints() const
Number of points supporting patch faces.
errorManip< error > abort(error &err)
A templated (N x N) square matrix of objects of <Type>, containing N*N elements, derived from Matrix.
Vector2D< Cmpt > x() const
Extract vector for row 0.
VectorSpace< Form, Cmpt, Ncmpts > normalised(const VectorSpace< Form, Cmpt, Ncmpts > &vs)
Representation of a 3D Cartesian coordinate system as a Vector of row vectors.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
A triangular face using a FixedList of labels corresponding to mesh vertices.
Vector2D< Cmpt > y() const
Extract vector for row 1.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
A 1D vector of objects of type <T>, where the size of the vector is known and can be used for subscri...
const dimensionedScalar e
Elementary charge.
static tmp< T > New(Args &&... args)
Construct tmp of T with forwarding arguments.
static const SymmTensor2D< Cmpt > zero
tensor rotationTensor(const vector &n1, const vector &n2)
Rotational transformation tensor from vector n1 to n2.
const Map< label > & meshPointMap() const
Mesh point map.
const labelList & meshPoints() const
Return labelList of mesh points in patch.
const word & constant() const
Return constant name.
Tensor< scalar > tensor
Tensor of scalars, i.e. Tensor<scalar>.
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]
bool set(const direction d) const
Is the vector in the direction d set.