CFD Online Discussion Forums

CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM Programming & Development (http://www.cfd-online.com/Forums/openfoam-programming-development/)
-   -   I know the cell ID, how do I find the cell coordinates? (http://www.cfd-online.com/Forums/openfoam-programming-development/110751-i-know-cell-id-how-do-i-find-cell-coordinates.html)

simpomann December 19, 2012 14:43

I know the cell ID, how do I find the cell coordinates?
 
Hey,

I have a complex case that seems to be crashing because of a single cell. Unluckily I cant identify it with the sets generated by checkMesh (allGeo allTopo). But I checked values immediately before the crash, so I think I know its ID (its nr 94495 in the list).

How do I find out where it is located? E.g. I need to find it's cell centre point or the vertices defining the cell (their coordinates).

So I know where inside my mesh the crash magic takes place.
Isn't there a cell list similiar to the face list in 0/polyMesh ?

My intention is to change meshing parameters in a way that this cell becomes normal.

Any help is highly appreciated! These crashes drive me absolutely crazy (as they happen after 5-20 hours of computation time and then its mesh-problem-hunt :eek:)

Greetings and thanks in advance,

Simon

Not sure if it is helpful but here is my checkMesh output (it's about tank filling).

Quote:

/*---------------------------------------------------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 2.1.0 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
Build : 2.1.0-bd7367f93311
Exec : checkMesh -allGeometry -allTopology
Date : Dec 19 2012
Time : 19:53:06
Host : "m-rocks.rz.fh-ingolstadt.de"
PID : 11425
Case : /home/ft8517/fin_2_191212
nProcs : 1
sigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).
fileModificationChecking : Monitoring run-time modified files using timeStampMaster
allowSystemOperations : Disallowing user-supplied system call operations
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time
Create polyMesh for time = 0
Time = 0
Mesh stats
points: 500821
faces: 1306654
internal faces: 1217023
cells: 409545
boundary patches: 16
point zones: 0
face zones: 0
cell zones: 0
Overall number of cells of each type:
hexahedra: 320740
prisms: 14208
wedges: 0
pyramids: 0
tet wedges: 1973
tetrahedra: 121
polyhedra: 72503
Checking topology...
Boundary definition OK.
Cell to face addressing OK.
Point usage OK.
Upper triangular ordering OK.
Face vertices OK.
Topological cell zip-up check OK.
<<Number of faces with non-consecutive shared points: 4. This might indicate a problem.
<<Writing 8 faces with non-standard edge connectivity to set edgeFaces
Number of regions: 1 (OK).
Checking patch topology for multiply connected surfaces ...
Patch Faces Points Surface topology Bounding box
walls 1 4 ok (non-closed singly connected) (2.08806 0.109125 0.0441667) (2.09158 0.111406 0.0441667)
inlet 44 55 ok (non-closed singly connected) (2.69246 0.897475 0.591929) (2.70942 0.900057 0.607066)
tank_patch0 21146 25372 ok (non-closed singly connected) (1.76868 -0.438435 -0.146001) (2.79999 0.633515 0.222356)
efr_patch0 16047 19241 ok (non-closed singly connected) (2.42331 0.42386 0.0460303) (2.71407 0.790794 0.585647)
entl_patch0 11192 13260 ok (non-closed singly connected) (2.48122 0.1231 0.0749478) (2.67755 0.721033 0.5623)
pistole_vorn_patch0 424 591 ok (non-closed singly connected) (2.66657 0.71442 0.51485) (2.68961 0.754364 0.558485)
pistole_hinten_patch01691 2104 ok (non-closed singly connected) (2.67042 0.75118 0.531703) (2.70942 0.900057 0.607466)
innenteile0_CATIASTL18319 20582 ok (non-closed singly connected) (1.86566 -0.364421 -0.1059) (2.77162 0.59733 0.2223)
innenteile1_CATIASTL7589 8529 ok (non-closed singly connected) (1.94061 0.270498 -0.146) (2.0825 0.491339 0.0681959)
innenteile2_CATIASTL3737 4198 ok (non-closed singly connected) (2.40454 -0.2589 0.0209999) (2.5714 -0.0595497 0.162726)
innenteile3_CATIASTL315 417 ok (non-closed singly connected) (2.49502 -0.269535 0.0209999) (2.53566 -0.141457 0.14723)
klappe_patch0 2633 3300 ok (non-closed singly connected) (2.39798 0.34564 0.0410414) (2.45511 0.435473 0.0850982)
flv_patch0 2892 3116 ok (non-closed singly connected) (2.53414 0.112304 0.12) (2.57011 0.148064 0.20453)
verschlusskopf_patch00 0 ok (empty)
einsatz_efr_patch0 3558 4600 ok (non-closed singly connected) (2.65035 0.696179 0.49897) (2.71363 0.790794 0.585647)
atmosphaere_patch0 43 78 ok (non-closed singly connected) (2.65798 0.772538 0.57277) (2.70486 0.782336 0.581323)
Checking geometry...
Overall domain bounding box (1.76868 -0.438435 -0.146001) (2.79999 0.900057 0.607466)
Mesh (non-empty, non-wedge) directions (1 1 1)
Mesh (non-empty) directions (1 1 1)
Boundary openness (-2.26301e-16 -1.04602e-16 2.10553e-15) OK.
Max cell openness = 1.50596e-15 OK.
Max aspect ratio = 47.3176 OK.
Minumum face area = 6.3327e-08. Maximum face area = 0.0012456. Face area magnitudes OK.
Min volume = 9.11448e-11. Max volume = 4.29548e-05. Total volume = 0.082352. Cell volumes OK.
Mesh non-orthogonality Max: 64.8266 average: 10.9586
Non-orthogonality check OK.
Face pyramids OK.
***Max skewness = 5.62295, 2 highly skew faces detected which may impair the quality of the results
<<Writing 2 skew faces to set skewFaces
Coupled point location match (average 0) OK.
***Error in face tets: 3 faces with low quality or negative volume decomposition tets.
<<Writing 3 faces with low quality or negative volume decomposition tets to set lowQualityTetFaces
Min/max edge length = 4.03237e-05 0.0367576 OK.
*There are 1591 faces with concave angles between consecutive edges. Max concave angle = 77.1632 degrees.
<<Writing 1591 faces with concave angles to set concaveFaces
Face flatness (1 = flat, 0 = butterfly) : average = 0.998086 min = 0.547233
*There are 142 faces with ratio between projected and actual area < 0.8
Minimum ratio (minimum flatness, maximum warpage) = 0.547233
<<Writing 142 warped faces to set warpedFaces
Cell determinant (wellposedness) : minimum: 5.70394e-05 average: 12.0743
***Cells with small determinant found, number of cells: 313
<<Writing 313 under-determined cells to set underdeterminedCells
***Concave cells (using face planes) found, number of cells: 22119
<<Writing 22119 concave cells to set concaveCells
Failed 4 mesh checks.
End

mturcios777 December 19, 2012 15:49

For cell center, you should be able to do mesh[cellI].x() to get the x-coordinate of the cell center. You can do a similar operation for y and z.

mulfal December 28, 2012 16:29

Hi Marco,
Good day and happy new year!
I have been trying to find the maximum value of the coordinates, say (max(y)) in order to calculate a hydrostatic pressure as:

hydrostatic pressure = rho * (g & (max(y) - mesh.C())),

where rho is the density of water and g is the gravitational acceleration.

I am having a problem, either in getting the maximum value of "y" from my computational domain ( which is the depth of the water level in principle, b/c I have air and water body in the domain) or the equation I am using is wrong. I have been trying to modify the postprocessing uttility (ptot), which calculates the total pressure (static + dynamic) in order to include the hydrostatic pressure but is giving me errors during compiliation. I have used a constant number say max(y) = 1, and still the solver is not compiling.

I have been thinking of finding max (y) when the density of the fluid is 1000 in the cell, which is the water density, in this way, I will be able to track the free surface dynamics and its corresponding height (max (y)). But this seems to be very difficult for the time being, for me at least. So I was thinking just to do in the simplest way by using constant value of the water depth though I couldn't do this one either.

I would appreciate if you can lend me a hand on how to access the maximum value of the coordinates and if the minus sign used in the above equation is correct in the programming point of view.

Thank you for your time and help.

With regards,

Mulualem

mturcios777 December 28, 2012 18:08

I think your problem is that you are mixing up types and operators. You are trying to subtract a single constant scalar value (max(y)) from a vector field (mesh.C()). I think you can agree this doesn't make sense. You may have to do the calculation cell by cell to extract the y value. Either change max(y) to a vector (create a temporary one that is just (0,1,0)*max(y)), and the subtraction and then dot product should work. As long as g is oriented in the y direction this should give you the answer you want.

mulfal December 30, 2012 05:53

Thank you Marco for your prompt reply.
You rightly pointed out my problem, it appears that it was mixing of types and operators. Now, the solver is working great with a temporary vector of the form "Foam::vector(0,1,0)". I will try the cell by cell method b/c that could help me to figure out the free surface dynamics by searching the cells with the maximum density.

Regards,

Mulualem

openfoammaofnepo February 25, 2014 13:11

Is it possible to extract the coordinates of the vertices belonging to one face in openbfoam? Assume the face is triangular, then the coordinates of the three vertices? Thank you.

Yeru June 30, 2014 10:29

Accessing coordinates in BC face by face
 
Dear all,

your advice is pretty pretty close to what I needed recently.
I'm trying to assign a scalar value to my boundaryField depending on each face's x- and z-value.

So since I'm doing it face by face in my 0/T-file figuratively I'm being handed a cell, I don't care about it's index, I just want to access THIS face's coordinates.

Is there something similar to Xcoord=mesh.x() I can use?

Here's my code with some placeholders
Code:

// 0/T-file
focus_inlet
{
 type          groovyBC;
variables    "circ=<value>;d_2=<value>;n=5;Xcoord=mesh.x();Zcoord=mesh.z();";
valueExpression "(<f(Xcoord,Zcoord)>    ?    1    :    0.1";
value        uniform 0.1;

EDIT: SOLVED! It's fairly simple
Code:

variables    "circ=<value>;d_2=<value>;n=5;Xcoord=pos().x;Zcoord=pos().z;";


All times are GMT -4. The time now is 00:54.