CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Post-Processing

Write cells and data intersecting a plane cuttingPlane

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree2Likes

Reply
 
LinkBack Thread Tools Display Modes
Old   April 8, 2008, 10:26
Default Hi Dragos Thank you very mu
  #41
ngj
Senior Member
 
Niels Gjoel Jacobsen
Join Date: Mar 2009
Location: Deltares, Delft, The Netherlands
Posts: 1,619
Rep Power: 25
ngj will become famous soon enoughngj will become famous soon enough
Hi Dragos

Thank you very much for sharing this tool ... I will definitely save a lot of disk space due to this way of exporting data (and a lot of time, since it would have taken me a long time to but this stuff together).
Though I have a minor question. Say that I am making a plane parallel to the x-y axes and I only want to export cells where x>3 and y<1, then do you know any clever way of removing the unwanted cells from the labelList?

Thanks

Niels
__________________
Please note that I do not use the Friend-feature, so do not be offended, if I do not accept a request.
ngj is offline   Reply With Quote

Old   April 9, 2008, 01:55
Default Hi Niels, I'm glad that you f
  #42
Senior Member
 
dmoroian's Avatar
 
Dragos
Join Date: Mar 2009
Posts: 647
Rep Power: 11
dmoroian is on a distinguished road
Hi Niels,
I'm glad that you found the code interesting.
I don't know the full answer to your question but this is a small guide line:
const labelList& cutCells = cutPlane.cells();
forAll(cutCells,cellI)
{
<blockquote>if(mesh.C()[cellI].x() <=3||mesh.c()[celli].y()>= 1)
<blockquote>//remove the item from the list
</blockquote>}</blockquote>

The problem is that I don't know how to remove one item from the list...
I've checked the class List and it doesn't seem to have functions for adding or removing items. Anyone knowing how to do it, please share your information.

Dragos
dmoroian is offline   Reply With Quote

Old   April 9, 2008, 04:16
Default Hi Dragos Yes, I thought of
  #43
ngj
Senior Member
 
Niels Gjoel Jacobsen
Join Date: Mar 2009
Location: Deltares, Delft, The Netherlands
Posts: 1,619
Rep Power: 25
ngj will become famous soon enoughngj will become famous soon enough
Hi Dragos

Yes, I thought of something similar to that approach, but came to the same question: How to remove items from the list.

Thanks for your time

- Niels
__________________
Please note that I do not use the Friend-feature, so do not be offended, if I do not accept a request.
ngj is offline   Reply With Quote

Old   April 9, 2008, 05:01
Default Hi Dragos, thanks, I will t
  #44
Senior Member
 
Fabian Braennstroem
Join Date: Mar 2009
Posts: 407
Rep Power: 10
braennstroem is on a distinguished road
Hi Dragos,

thanks, I will try the vtk write as soon as I know how to do it ;-)

Maybe, you could help in the meanwhile with a different problem. As soon as I have the subsette/plane, is there a chance to find the maximum or minimum of a variable and access the cell id and the absolute coordinates?

Fabian
braennstroem is offline   Reply With Quote

Old   April 9, 2008, 12:08
Default Hi Dragos, thanks, though i
  #45
Senior Member
 
Fabian Braennstroem
Join Date: Mar 2009
Posts: 407
Rep Power: 10
braennstroem is on a distinguished road
Hi Dragos,

thanks, though it seems to be a global max or min. One probably has to narrow the variable values to the subsetter mesh!?

Fabian
braennstroem is offline   Reply With Quote

Old   April 9, 2008, 15:44
Default Ok, try: max(scalarFlds)
  #46
Senior Member
 
dmoroian's Avatar
 
Dragos
Join Date: Mar 2009
Posts: 647
Rep Power: 11
dmoroian is on a distinguished road
Ok, try: max(scalarFlds[i])
dmoroian is offline   Reply With Quote

Old   April 10, 2008, 07:53
Default Thanks, I just try something l
  #47
Senior Member
 
Fabian Braennstroem
Join Date: Mar 2009
Posts: 407
Rep Power: 10
braennstroem is on a distinguished road
Thanks, I just try something like:

label buffmin = findMin(scalarFlds1);
label buffmax = findMax(scalarFlds1);
forAll(scalarFlds1, i)
{
buffmin = findMin(scalarFlds1[i]);
buffmax = findMax(scalarFlds1[i]);
}

but I am not sure yet, if really the selected plane is checked... right now, it does not look like it; I assume that it gets the min and max of the last iteration in the forall loop.

Fabian
braennstroem is offline   Reply With Quote

Old   April 10, 2008, 08:17
Default Hi Fabian, Of course your loo
  #48
Senior Member
 
dmoroian's Avatar
 
Dragos
Join Date: Mar 2009
Posts: 647
Rep Power: 11
dmoroian is on a distinguished road
Hi Fabian,
Of course your loop will give you only the last list maximum/minimum. In my opinion it should be done like this:
labelList buffmin(scalarFlds.size());
labelList buffmax(scalarFlds.size());
forAll(scalarFlds,i)
{
<blockquote>buffmin[i] = findMin(scalarFlds[i]);
buffmax[i] = findMax(scalarFlds[i]);
}</blockquote>

or instead of the loop (I'm not sure if it works):
buffmin = findMin(scalarFlds);
buffmax = findMax(scalarFlds)

But you definitely need a list for both minimum label and maximum label.

Dragos
dmoroian is offline   Reply With Quote

Old   April 10, 2008, 13:45
Default Thanks, I tried it and this ap
  #49
Senior Member
 
Fabian Braennstroem
Join Date: Mar 2009
Posts: 407
Rep Power: 10
braennstroem is on a distinguished road
Thanks, I tried it and this approach:

label buffmin = findMin(scalarFlds1);
label buffmax = findMax(scalarFlds1);

//labelList buffmin(scalarFlds1.size());
//labelList buffmax(scalarFlds1.size());
label max = 0.0;
label min = 0.0;
forAll(scalarFlds1, i)
{
buffmin = findMin(scalarFlds1[i].internalField());
buffmax = findMax(scalarFlds1[i].internalField());

if (p[buffmin] < min)
{
min=p[buffmin];
}
if (p[buffmax] > max)
{
max=p[buffmax];
}
}

unfortunately with no success. Maybe there is a problem with the element index. 'buffmin' is the local index of the local min pressure cell, but p[buffmin] gives the pressure at the global index, which is different to the local index?! Am I still running in the wrong direction?

Fabian
braennstroem is offline   Reply With Quote

Old   April 10, 2008, 16:43
Default First of all, it seems that I
  #50
Senior Member
 
dmoroian's Avatar
 
Dragos
Join Date: Mar 2009
Posts: 647
Rep Power: 11
dmoroian is on a distinguished road
First of all, it seems that I don't understand you very well. If buffmin is of type label consequently scalarFlds1 must have only one element so there is no point in the forAll loop.
Second, it is dangerous to set the minimum min = 0 as well the maximum max = 0. What if you have all relative pressures larger than 0, or what if all of them are less than 0. In either case you will end up with a wrong min or max.
Anyway if you have only one scalar in scalarFlds1, let say pressure p, then the index of the min pressure should be given by
min = findMin(scalarFlds1); pmin = scalarFlds1[min];
and similarly the max pressure sould be given by
max = findMax(scalarFlds1); pmax = scalarFlds1[max];
dmoroian is offline   Reply With Quote

Old   April 11, 2008, 05:15
Default Hi Dragos, it works :-) Act
  #51
Senior Member
 
Fabian Braennstroem
Join Date: Mar 2009
Posts: 407
Rep Power: 10
braennstroem is on a distinguished road
Hi Dragos,

it works :-) Actually min and max = 0 were just a 'dummy'.
I use these lines with gMin, gMax (as in bound.c):

label buffmin = findMin(scalarFlds1);
label buffmax = findMax(scalarFlds1);
label fmin = findMin(scalarFlds1[buffmin].internalField());
label fmax = findMax(scalarFlds1[buffmin].internalField());

Info << "size of scalarFlds: " << scalarFlds1.size() << " " << endl; // nur 1
Info << "size of scalarFlds: " << scalarFlds1[buffmin].internalField().size() << " " << endl; // nur 1

//Info << "local index of min: " << buffmin << " local index of max: " << buffmax << " " << endl; //0
Info << "local index of min: " << fmin << " local index of max: " << fmax << " " << endl;

Info << "min: " << gMin(scalarFlds1[buffmin].internalField()) << endl;
Info << "max: " << gMax(scalarFlds1[buffmin].internalField()) << endl;
Info << "average: " << gAverage(scalarFlds1[buffmin].internalField()) << endl;
//Info << "coords of min in global mesh: " << mesh.C()[fmin] << endl;
Info << "coords of min in local mesh: " << subsetter1.subMesh().C()[fmin] << endl;
//Info << "coords of max in global mesh: " << mesh.C()[fmax] << endl;
Info << "coords of max in local mesh: " << subsetter1.subMesh().C()[fmax] << endl;


Once again, thanks a lot for your help!
Fabian
braennstroem is offline   Reply With Quote

Old   May 21, 2008, 08:48
Default Hi all I have been using cu
  #52
ngj
Senior Member
 
Niels Gjoel Jacobsen
Join Date: Mar 2009
Location: Deltares, Delft, The Netherlands
Posts: 1,619
Rep Power: 25
ngj will become famous soon enoughngj will become famous soon enough
Hi all

I have been using cuttingPlane as the one give by Dragos with a lot of success, but unfortunately the there is a snake in Paradise... When I write the cuttingPlane mesh using the command:

subsetter.subMesh().write();

it is written to a polyMesh-directory in the current time-directory. This yields three problems:

1. When I want to view the datasets of the entire domain, which are written less often, ParaFoam uses the latest polyMesh directory to generate the mesh, and since there are fewer points in this mesh than in the complete mesh, the data and mesh are incompatible and thus paraFoam terminates.

2. A related problem occur when using latestTime as startTime. Apparently all time directories are searched and the latest polyMesh file is used, leading again to a mismatch in the number of cells in the polyMesh and the field data, e.g. U.

3. I am also going to due moving meshes at some point. There I need to output the entire mesh after the updating, and I can imagine, that one of the polyMesh-directories is then overwritten by the other.

Thus my question is, if there is any way to set a sub-directory to the current time-directory in which all cutting-plane data is written? Looking into the fvMesh in Doxygen does not seem to allow such. The crude solution is to make a system call and generate a folder and move all cuttingPlane data at runTime, but it is very inelegant.

Thanks a lot for your help,

Niels
__________________
Please note that I do not use the Friend-feature, so do not be offended, if I do not accept a request.
ngj is offline   Reply With Quote

Old   May 22, 2008, 07:08
Default Hello Niels, You are right, t
  #53
Senior Member
 
dmoroian's Avatar
 
Dragos
Join Date: Mar 2009
Posts: 647
Rep Power: 11
dmoroian is on a distinguished road
Hello Niels,
You are right, there are small issues regarding the mesh.
I did not have enough time to look at it, so the fast workaround was to reconstruct an new case directory as a postprocessing phase. This new directory contains only the interesting data (either 3D or the cutting plane). The main steps:
# 1. generate the new case directory

# 2. add BC in processorX/constant/boundary for every other processor at a time

# 3. remove the probes from each controlDict file

# 4. merge meshes

# 5 map fields

I hope this is useful,
Dragos
dmoroian is offline   Reply With Quote

Old   May 26, 2008, 06:55
Default Hi Dragos That might be wha
  #54
ngj
Senior Member
 
Niels Gjoel Jacobsen
Join Date: Mar 2009
Location: Deltares, Delft, The Netherlands
Posts: 1,619
Rep Power: 25
ngj will become famous soon enoughngj will become famous soon enough
Hi Dragos

That might be what I am looking for. I will take a look at it, and then I might return to you.

Thanks for the help.

- Niels
__________________
Please note that I do not use the Friend-feature, so do not be offended, if I do not accept a request.
ngj is offline   Reply With Quote

Old   June 3, 2008, 12:29
Default Hi, I used the following co
  #55
Senior Member
 
John Deas
Join Date: Mar 2009
Posts: 160
Rep Power: 8
johndeas is on a distinguished road
Hi,

I used the following code to define a subset:

point pnt(6.251,0.501,0.001);
vector direction(1,0,0);
plane pl1(pnt,direction);
cuttingPlane cutPlane(mesh,pl1);
const labelList& cutCells = cutPlane.cells();
Info << cutCells << endl;

when I launch my calculation on one processor I retrieve a list of intersecting cells, but when I launch it in parallel I get 0(). Has someone already faced this, using cuttingPlane in parallel ?
johndeas is offline   Reply With Quote

Old   June 4, 2008, 02:38
Default John, If you replace 'Info'
  #56
Senior Member
 
Mark Olesen
Join Date: Mar 2009
Location: http://olesenm.github.io/
Posts: 777
Rep Power: 18
olesen will become famous soon enough
John,

If you replace 'Info' with 'Pout' you should see that cuttingPlane does work, but just didn't happen to cut any cells on the first cpu.

If you need to recombine the parallel lists of cells, you'll need to use something like ListListOps::combineOffset on the master process.
I don't know offhand where to look, but grep through some of the source to find something equivalent.
olesen is offline   Reply With Quote

Old   June 4, 2008, 02:39
Default hi John I have not experien
  #57
ngj
Senior Member
 
Niels Gjoel Jacobsen
Join Date: Mar 2009
Location: Deltares, Delft, The Netherlands
Posts: 1,619
Rep Power: 25
ngj will become famous soon enoughngj will become famous soon enough
hi John

I have not experienced such problems, and I have adopted the exact way described by Dragos above. If you are interested, I can send you my files, so you can try them.

Best regards,

Niels
__________________
Please note that I do not use the Friend-feature, so do not be offended, if I do not accept a request.
ngj is offline   Reply With Quote

Old   July 22, 2008, 10:13
Default Hi Dragos I have begun look
  #58
ngj
Senior Member
 
Niels Gjoel Jacobsen
Join Date: Mar 2009
Location: Deltares, Delft, The Netherlands
Posts: 1,619
Rep Power: 25
ngj will become famous soon enoughngj will become famous soon enough
Hi Dragos

I have begun looking into the merging, and I have tried to follow the link above.
Though, it is still somewhat not obvious how to merge the meshes generated by, say 4 processors.
I hope you can find time to elaborate a little on the procedure.
I have a significant amount of time steps where I need to merge and map the data, thus I am thinking of making something generic, which I of course will be happy to share, when it is done.

Thanks for all your help so far,

Niels
__________________
Please note that I do not use the Friend-feature, so do not be offended, if I do not accept a request.
ngj is offline   Reply With Quote

Old   August 11, 2008, 02:54
Default Hy Niels, Sorry for the delay
  #59
Senior Member
 
dmoroian's Avatar
 
Dragos
Join Date: Mar 2009
Posts: 647
Rep Power: 11
dmoroian is on a distinguished road
Hy Niels,
Sorry for the delay... I had a long and nice vacation, but the internet connection was very bad.
Well, indeed it is tricky to merge a decomposed mesh with mergeMeshes. There are several answers on this forum regarding only this problem. The sensitive part is that you need to modify the boundary file in each processor directory.
Compared to producing the slice, putting back together the bits from a parallel computation is a lot tougher. The solution I adopted is ugly although is general, and comprises both perl and bash scripts.
Due to some secrecy conditions that I have to comply with, I cannot publicly share this files.
But there is absolutely no problem to talk or share the information either here or on private, and then come up with a tool on wiki.

Regards,
Dragos
dmoroian is offline   Reply With Quote

Old   August 21, 2008, 08:34
Default Hi Dragos and any other who is
  #60
ngj
Senior Member
 
Niels Gjoel Jacobsen
Join Date: Mar 2009
Location: Deltares, Delft, The Netherlands
Posts: 1,619
Rep Power: 25
ngj will become famous soon enoughngj will become famous soon enough
Hi Dragos and any other who is interested

I have just returned from vacation and I have begun working on the merging/mapping. I have made a small shell script which takes care of the merging, and I believe it to be generic.
For all you others, then I can tell that I have taken direct contact with Dragos, so hopefully there will soon be a nicely wrapped up tool ready to public use.

Enjoy the rest of the day (and the day after as well...),

Niels
__________________
Please note that I do not use the Friend-feature, so do not be offended, if I do not accept a request.
ngj is offline   Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
UDF to write a data file lichun Dong FLUENT 2 August 4, 2005 11:21
write transient data of plane Fabian CFX 2 June 27, 2005 02:38
UDF TO WRITE OUT DATA!! Aishwarya FLUENT 1 August 14, 2004 10:32
UDF TO WRITE OUT DATA!! Aishwarya FLUENT 0 August 12, 2004 10:11
UDF to read and write data Mcgregor FLUENT 4 June 9, 2003 13:21


All times are GMT -4. The time now is 16:15.