
[Sponsors] 
June 2, 2010, 05:03 
Convert coordinate system in OpenFOAM

#1 
Senior Member
Gijsbert Wierink
Join Date: Mar 2009
Posts: 383
Rep Power: 10 
Dear all,
I am working with a project where we would like to implement a surface capturing method that requires information of the coordinate system in the form x(i,j,k). I have understood that OpenFOAM uses a cell list of the form x(i), so you do not really know the ijk position of a cell. The grid we want to use is a simple orthogonal grid. I would like to know the following:
__________________
Regards, Gijs 

June 2, 2010, 05:31 

#2 
Senior Member
Niels Gjoel Jacobsen
Join Date: Mar 2009
Location: Deltares, Delft, The Netherlands
Posts: 1,676
Rep Power: 26 
Hi
Well, you can achieve the information from the following building blocks: Code:
// Neighbour and owner cells to given face const unallocLabelList & nei = mesh.neighbour(); const unallocLabelList & own = mesh.owner(); forAll( mesh.C(), celli) { // This is a list of faces of which the cell is build const cell& cProp = mesh_.cells()[celli]; // Loop over all faces forAll(cProp, facei) { // Using own and nei you can determine x(i,j,k) } } Best regards, Niels 

June 2, 2010, 05:45 

#3 
Senior Member
Gijsbert Wierink
Join Date: Mar 2009
Posts: 383
Rep Power: 10 
Hi Niels,
Fantastic! Thank you for your quick reply! The method I want to implement is called the Constrained Interpolation Profile (CIP) method. It is basically a subgrid reconstruction method which you can use to capture an interface. I am interested in modelling interaction between gasliquid interfaces and particles (or drops of a third phase/species). By using CIP on an orthogonal grid I hope to be able to implement more coupled physics easier, like Marangoni effects, EhpH, electrophoretic forces etc. Apologies for the noobquestion , but would you have some hints on how to determine x(i,j,k) from own and nei? Thanks in advance!
__________________
Regards, Gijs 

June 2, 2010, 05:52 

#4 
Senior Member
Niels Gjoel Jacobsen
Join Date: Mar 2009
Location: Deltares, Delft, The Netherlands
Posts: 1,676
Rep Power: 26 
Hi Gijs
Okay, that is interesting. I did an implementation of a DDR scheme following Harvie & Fletcher earlier this year. It has some better properties with respect to keeper a clean interface, however my implementation unfortunately seems to be somewhat faulty. With respect to neighbouring/owner, then Code:
// Cell labels label owner = own[cProp[facei]]; label neigh = nei[cProp[facei]]; if (owner == celli) // then the neighbour is actaully neigh else // neigh == celli // then the neighbour is actually owner Best regards, Niels 

June 3, 2010, 07:44 
How to check alpha1 value on each side of a face and "grow"

#5 
Senior Member
Gijsbert Wierink
Join Date: Mar 2009
Posts: 383
Rep Power: 10 
Hi Niels,
Thanks for your answer. Due to some other issues (related to the physicochemistry we want to solve) we decided to use only part of the CIP method, so we don't really need the x(i,j,k) coordinates at this point. Sorry for the hassle about that . To implement the (partial) CIP method I want to try transporting a solid using interFoam first, as a test. I intend to kind of freeze the solid and liwuid around it and then translate it according to the flow field. This is to prevent diffussion of the solid (which will otherwise happen with this VOF method). In order to do this I need to browse through all the faces and check whether there is on one cell center water and on the other side air. If I then find the faces that are "interface faces" I want to "grow" a (few) cell layer(s) around the solid object and then transport the "grown solid block". Then, I want to integrate the forces on the faces of the grown interface. Would you have some hints how to implement this? I have some pseudocode below, but lots is still missing ... Code:
// Neighbour and owner cells to given face const unallocLabelList& nei = mesh.neighbour(); const unallocLabelList& own = mesh.owner(); // Loop through all cell faces forAll(mesh.Sf(), facei) { // Check for all faces whether it has on one side // a cell center with alpha1 = 0 and a cell center // with alpha1 = 1 on the other side. forAll( mesh.C(), celli) { // Grow a cell layer in the direction of alpha1 = 0 } }
__________________
Regards, Gijs 

June 7, 2010, 07:58 

#6 
Senior Member
Niels Gjoel Jacobsen
Join Date: Mar 2009
Location: Deltares, Delft, The Netherlands
Posts: 1,676
Rep Power: 26 
Hi Gijs
Unfortunately I cannot help you with your problem, however if the object you are going to consider remain a simple geometry throughout the computation, then you might want to consider moving the boundary between the two fluids, e.g. [1]. Z. Tukovic has further developed a FAMroutine, where you solve given differential equations on a surface. As you have chemistry in your problem it might be relevant so forth some of the reactions are limited to take place on the interface. Good luck, Niels [1]: Tukovic, Z. and Jasak, H. (2007). Updated Lagrangian finite volume solver for large deformation dynamic response of elastic body. Transactions of FAMENA, 31(1), 55–70. 

June 7, 2010, 08:42 

#7  
Senior Member
Gijsbert Wierink
Join Date: Mar 2009
Posts: 383
Rep Power: 10 
Hi Niels,
Thanks for your help. Quote:
__________________
Regards, Gijs 

June 7, 2010, 08:51 

#8 
Senior Member
Niels Gjoel Jacobsen
Join Date: Mar 2009
Location: Deltares, Delft, The Netherlands
Posts: 1,676
Rep Power: 26 
You are welcome. I remembered having a link to the mentioned article, see [1]. Especially section 5 might interest you.
Further, regarding FAMframework it is already implemented in 1.5dev and there is a tutorial where the dispersion of a species on a curved surface is modelled. Best regards, Niels [1] http://powerlab.fsb.hr/ped/kturbo/Op...tionFAMENA.pdf 

June 17, 2010, 03:37 

#9 
Senior Member
Gijsbert Wierink
Join Date: Mar 2009
Posts: 383
Rep Power: 10 
Hi Niels,
Is there an easy way to loop through all (internal) faces in a mesh? Something like the following perhaps (although this doesn't compile)? Code:
const unallocLabelList& listOfFaces = mesh.faces(); forAll(listOfFaces, faceI) { // do stuff }
__________________
Regards, Gijs 

June 17, 2010, 06:02 

#10 
Senior Member
Niels Gjoel Jacobsen
Join Date: Mar 2009
Location: Deltares, Delft, The Netherlands
Posts: 1,676
Rep Power: 26 
Hi Gijs
I have looped over the internal faces of a surfaceScalarField, e.g. phi, when doing this. I guess you need to flux anyway for your two phase solver, so that might be a way to go. Best regards, Niels 

June 17, 2010, 09:01 

#11 
Senior Member
Gijsbert Wierink
Join Date: Mar 2009
Posts: 383
Rep Power: 10 
Hi Niels,
Right, so you're suggesting to interpolate the volScalarField onto the faces and then loop through them? I mean, phi is already the face flux, isn't it? I am dealing with cell centered values (alpha). I guess it can be done, but is it efficient?
__________________
Regards, Gijs 

June 17, 2010, 09:07 

#12 
Senior Member
Niels Gjoel Jacobsen
Join Date: Mar 2009
Location: Deltares, Delft, The Netherlands
Posts: 1,676
Rep Power: 26 
What I mean is that when I have needed what you are doing, I was considering the direction of the flux over the cell faces, hence I needed to consider the sign of phi on every face relative to the orientation of the hence. This made it a obvious choice.
Could you elaborate on what it is you need? Best regards, Niels 

June 17, 2010, 09:22 

#13 
Senior Member
Gijsbert Wierink
Join Date: Mar 2009
Posts: 383
Rep Power: 10 
Sure, sorry for the fuzzyness .
I have a phase field alpha1, which is set to 1 inside a small blob and 0 in the main fluid. What I want to do is "pad" the blob with another phase, which is actually the same as the main fluid, just name differently for convenience. Let's call that phase beta. So, I want to loop through all faces and compare alpha1 on each side. When alpha1 is significantly different on one side than the other, it is an interface and I want to "pad" the adjacent main fluid as beta. Now I have this: Code:
forAll(alpha1, cellI) { beta[cellI] = alpha1[cellI]; label ownerCell = mesh.faceOwner()[cellI]; label neighbourCell = mesh.faceNeighbour()[cellI]; // then, compare owner and neighbour }
__________________
Regards, Gijs 

June 17, 2010, 11:33 
Jep, face flux did the trick!

#14  
Senior Member
Gijsbert Wierink
Join Date: Mar 2009
Posts: 383
Rep Power: 10 
Hi Niels,
I got rid of my strawberry! Quote:
__________________
Regards, Gijs 

June 18, 2010, 08:16 

#15 
Senior Member
Niels Gjoel Jacobsen
Join Date: Mar 2009
Location: Deltares, Delft, The Netherlands
Posts: 1,676
Rep Power: 26 
You are welcome.
Have a good weekend, Niels 

November 19, 2010, 05:32 

#16 
Member
Erik Arlemark
Join Date: Mar 2009
Location: Eindhoven, Netherlands
Posts: 47
Rep Power: 9 
Hi all and Gijsber,
I would like to find the minimum and maximum x, y and z coordinates in my geometry while running my solver. I can see these values from the case's blockMeshDict file but I would like to avoid manually entering these scalar values by hand in the solver code. I find that when NOT running in parallel the following lines works fine: //values below somewhere in the middle of the domain scalar minX=0, minY=0, minZ=0, maxX=0, maxY=0, maxZ=0; forAll( mesh.cellCentres(), cellPrev) { const cell& cProp = mesh.cells()[cellPrev]; forAll(cProp, facei) { if (mesh.faceCentres()[cProp[facei]].x() < minX) minX=mesh.faceCentres()[cProp[facei]].x(); if (mesh.faceCentres()[cProp[facei]].y() < minY) minY=mesh.faceCentres()[cProp[facei]].y(); if (mesh.faceCentres()[cProp[facei]].z() < minZ) minZ=mesh.faceCentres()[cProp[facei]].z(); if (mesh.faceCentres()[cProp[facei]].x() > maxX) maxX=mesh.faceCentres()[cProp[facei]].x(); if (mesh.faceCentres()[cProp[facei]].y() > maxY) maxY=mesh.faceCentres()[cProp[facei]].y(); if (mesh.faceCentres()[cProp[facei]].z() > maxZ) maxZ=mesh.faceCentres()[cProp[facei]].z(); } } However, as I run these lines in a parallel simulation they do not yield the min and max coordinate values for the respective processor domain. I take it that this is because some processors do not own their boundary face, but a cell of the neighboring processor cell does, and therefor these bordering vertices are not seen. Any suggestion for how to retrieve the min and max coordinate values is appreciated! Cheers, Erik Arlemark 

November 19, 2010, 06:09 

#17 
Senior Member
Laurence R. McGlashan
Join Date: Mar 2009
Posts: 370
Rep Power: 15 
If the mesh doesn't change during your simulation then the easiest thing to do is to have a separate 'preprocessing' application that works things like this out for you in serial. Otherwise you have to do something more complicated.
__________________
Laurence R. McGlashan :: Website 

November 20, 2010, 05:20 

#18 
Member
Erik Arlemark
Join Date: Mar 2009
Location: Eindhoven, Netherlands
Posts: 47
Rep Power: 9 
Thank you Laurence R. McGlashan,
I chose to get the coordinates by creating an extra dictionary in my case, where I put in all the user inputs being specific for my case. Cheers, Erik 

September 12, 2011, 10:08 

#19 
Member
Santiago
Join Date: Dec 2009
Posts: 85
Rep Power: 8 
Hi all.
I'm working in developing a fluid structure interaction solver for solidfluid. I was looking into obtaining the coordinates of each of the cells at the interface so this post helps. While searching the forums I notice how people seem to know all kind of classes in openfoam and what they actually do. My big question is: Is there a list or guide or something where I can learn about the different classes that I can use in OF? At this point I would like to know about class that operate on the mesh. I will appreciate any help. Thanks 

Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
spherical coordinate system in Fluent  Amir  FLUENT  1  March 8, 2015 13:25 
Defining local coordinate system  bastil  OpenFOAM PreProcessing  11  August 8, 2011 08:26 
Modified OpenFOAM Forum Structure and New MailingList  pete  Site News & Announcements  0  June 29, 2009 05:56 
Adventure of fisrst openfoam installation on Ubuntu 710  jussi  OpenFOAM Installation  0  April 24, 2008 14:25 
OpenFOAM Debian packaging current status problems and TODOs  oseen  OpenFOAM Installation  9  August 26, 2007 13:50 