[virginie_e]

Hi Foamers,

I have a question about OpenFOAM and the mesh handling.

I am trying to code an interTrackFoam solver which would keep the quality of the mesh above a certain level (typically, I use as an indicator the ratio of the inradius and the circumradius of the tetrahedra of the mesh and I do not allow it to be under a minimum value).

So I coded functions which enables the cells of the mesh to check this criterion and when I look at the resulting mesh, everything seems fine: there is no irregularity or holes in the mesh and the process runs without any error message.

However, when I perform this function, the calculus of U and p goes on a little while but quickly produce a time step continuity error.

I think this is due to the fact that I did not correctly orient my faces: I read somewhere that an internal face should be such has neighbour label > owner label and such that the ordering of the points produce a normal which points into the neighbour cell (with the right-hand rule).

Is there a way to set the orientation of the faces right automatically with the class directTopoChange? I saw variables such as flipFaceFlux... Is that it? Should I set this to true?

Or have you got any idea of what this error could be due to? Thank you a lot in advance.

Virginie

[deepsterblue]

I'm quite sure that the error is due to a wrong flux value, and you're right - flipFaceFlux is used to reverse the face-flux if the face changed orientation due to a topo-change. However, you can choose to correct fluxes explicitly using an additional pressure Poisson equation - take a cue from icoDyMFoam, where a correctPhi.H is used.

It's up to you to decide whether the Poisson solution is warranted - it may or may not be a performance issue. For simple topo-changes, flipFaceFlux is a convenient way to maintain divergence-free fluxes, but more complicated situations might simply be too difficult to handle, and a Poisson solution is a more convenient solution (which I'm really not fond of).

Hope this helps.

[virginie_e]

Thank you very much Sandeep for your answer. It helped me indeed to solve the problems I had. I had though another question: I read somewhere that for a given face, the cell label of the neighbour cell should be greater than the label owner. What might happen if it is not the case? If i set the option cells renumbering true in the directTopoChange function update mesh, does it set cell labels so that this rule is respected?

Thank you a lot.

Virginie

[deepsterblue]

Yup. It's good to be consistent about owner-neighbour ordering, since this approach has the advantage of being able to assume (without a priori knowledge about the mesh) that face normals will point from owner to neighbour, and that face right-handedness defines that direction as well. The face-normal direction is fairly important while defining gradient/divergence operators, and I suspect that improper owner-neighbour numbering won't have an impact on solution results, but I haven't tested this. I think it would be easier to play safe and stick to the convention.

A cursory glance at the directTopoChange code suggests that the owner-neighbour ordering is taken care of during the re-ordering stage. I found this bit of code in directTopoChange.C:

            // Renumber owner/neighbour. Take into account if neighbour suddenly
            // gets lower cell than owner.
            forAll(faceOwner_, faceI)
            {
                label own = faceOwner_[faceI];
                label nei = faceNeighbour_[faceI];

                if (own >= 0)
                {
                    // Update owner
                    faceOwner_[faceI] = localCellMap[own];

                    if (nei >= 0)
                    {
                        // Update neighbour.
                        faceNeighbour_[faceI] = localCellMap[nei];

                        // Check if face needs reversing.
                        if
                        (
                            faceNeighbour_[faceI] >= 0
                         && faceNeighbour_[faceI] < faceOwner_[faceI]
                        )
                        {
                            faces_[faceI] = faces_[faceI].reverseFace();
                            Swap(faceOwner_[faceI], faceNeighbour_[faceI]);
                        }
                    }
                }
                else if (nei >= 0)
                {
                    // Update neighbour.
                    faceNeighbour_[faceI] = localCellMap[nei];
                }
            }

Hope this helps.