[ngj]

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

[dmoroian]

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

[ngj]

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

[braennstroem]

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]

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

[dmoroian]

Ok, try: max(scalarFlds[i])

[braennstroem]

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

[dmoroian]

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

[braennstroem]

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

[dmoroian]

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];

[braennstroem]

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

[ngj]

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

[dmoroian]

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

[ngj]

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

[johndeas]

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 ?

[olesen]

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.

[ngj]

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

[ngj]

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

[dmoroian]

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

[ngj]

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