CFD Online Discussion Forums - Wrong behaviour when appending vector to vectorField?

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)

- OpenFOAM Bugs (https://www.cfd-online.com/Forums/openfoam-bugs/)

- - Wrong behaviour when appending vector to vectorField? (https://www.cfd-online.com/Forums/openfoam-bugs/153748-wrong-behaviour-when-appending-vector-vectorfield.html)

Wrong behaviour when appending vector to vectorField?

This code:

vectorField test;
vector t1(1.0,2.0,3.0);
test.append(t1);
vector& t2 = test[0];
test.append(t2);
Info << "test = " << test << endl;

results in the following output:

test = 2((1 2 3) (0 2 3))

The first component of the appended vector is not appended correctly.

The error does not occur if t2 is a copy of test[0] rather than a reference to it.

Here is the situation where I discovered the bug:

//Calculate a pointField and if it only contains a single point
//append this, so it occurs twice in the pointField
pointField pf;
someFunctionPopulatingPointField(pf);
if (pf.size()==1)
{
pf.append(pf[0]);
}

The bug was reported here:
http://www.openfoam.org/mantisbt/view.php?id=1729#c4882
It was marked as "resolved" by noting that one should work with a copy of instead of a reference to test[0]. As far as I can see this is not a fix but a work-around.

Could some C++ expert comment on whether this is a genuine bug that should be properly fixed or if am I just doing dodgy stuff in my coding?

Hi,

I would call this "undefined behaviour". If you take a look at how append method is implemented (Field is a child of List):

Code:

template<class T>

inline void Foam::List<T>::append(const T& t)

{

    setSize(size()+1, t);

}

and

Code:

template<class T>

void Foam::List<T>::setSize(const label newSize, const T& a)

{

    label oldSize = label(this->size_);

    this->setSize(newSize);



    if (newSize > oldSize)

    {

        register label i = newSize - oldSize;

        register T* vv = &this->v_[newSize];

        while (i--) *--vv = a;

    }

}

So append first resizes list and then set last element to the new value.

Now let's look at setSize implementation (src/OpenFOAM/containers/Lists/List/List.C:318), below is just relevant part

Code:

template<class T>

void Foam::List<T>::setSize(const label newSize)

{

    ...

    if (newSize != this->size_)

    {

        if (newSize > 0)

        {

            T* nv = new T[label(newSize)];



            if (this->size_)

            {

                ...copy old elements into newly allocated memory...

            if (this->v_) delete[] this->v_;



            this->size_ = newSize;

            this->v_ = nv;

        }

        ...

    }

}

After setSize call your t2 vector points to invalid location (memory was freed). To fix this one can a) use copy instead of reference, b) check if list contains a (in setSize(const label newSize, const T& a)). Variant a is easier ;)

Thanks a lot for that useful reply, alexeym!

I was not aware that appending to vectorFields and pointFields copied the whole field into new memory.

The right thing to do in my particular case would probably be to work instead with a DynamicList<vector> (with an appropriately guessed initial size) instead of a vectorField. This way I'll avoid all that copying.

Hi,

Well, not quite sure you can avoid copying with DynamicList:

Code:

template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>

inline void Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::setSize

(

    const label nElem

)

{

        ...

        {

            capacity_ = max

            (

                nElem,

                label(SizeInc + capacity_ * SizeMult / SizeDiv)

            );

        }



        List<T>::setSize(capacity_);

    }

    

    // adjust addressed size

    List<T>::size(nElem);

}

Though since DynamicList grows in chunks, memory motion events happen less frequently.

Yes, hence my comment in brackets: "(with an appropriately guessed initial size)".

If for instance you know that your pointField/DynamicList<point> will be a subset of the points of a particular face in your mesh, just make sure to allocate memory for 'nFacePoints' points from the start to avoid copying.

Cheers,
Johan