invalid read in multi-phase

[cwsmith]

In solver.inp I have set:

 Solve Level Set : 2                   # =>iLSet
 Body Force Option: Vector
 Body Force: 0.0 0.0 0.0   

When the 'Body Force' vector is parsed here:
https://github.com/PHASTA/phasta/blame/06d45c3e23cf9519431c2d5c7bbad33ddf2fbe7a/phSolver/common/input_fform.cc#L484-L486
vec is accessed well beyond the size of vector of length 3. Valgrind reports this as an invalid read:

==21108== Invalid read of size 8
==21108==    at 0x46509E: input_fform(char*) (input_fform.cc:475)
==21108==    by 0x448698: phasta (phasta.cc:138)
==21108==    by 0x4483E7: main (main.c:7)
==21108==  Address 0x6dea3d8 is 0 bytes after a block of size 24 alloc'd
==21108==    at 0x4C2695A: operator new(unsigned long) (vg_replace_malloc.c:298)
==21108==    by 0x473A93: __gnu_cxx::new_allocator<double>::allocate(unsigned long, void const*) (new_allocator.h:104)
==21108==    by 0x4736B1: __gnu_cxx::__alloc_traits<std::allocator<double> >::allocate(std::allocator<double>&, unsigned long) (alloc_traits.h:182)
==21108==    by 0x4730F5: std::_Vector_base<double, std::allocator<double> >::_M_allocate(unsigned long) (stl_vector.h:170)
==21108==    by 0x472A25: double* std::vector<double, std::allocator<double> >::_M_allocate_and_copy<__gnu_cxx::__normal_iterator<double const*, std::vector<double, std::allocator<double> > > >(unsigned long, __gnu_cxx::__normal_iterator<double const*, std::vector<double, std::allocator<double> > >, __gnu_cxx::__normal_iterator<double const*, std::vector<double, std::allocator<double> > >) (stl_vector.h:1224)
==21108==    by 0x472176: std::vector<double, std::allocator<double> >::operator=(std::vector<double, std::allocator<double> > const&) (vector.tcc:195)
==21108==    by 0x462BA7: input_fform(char*) (input_fform.cc:312)
==21108==    by 0x448698: phasta (phasta.cc:138)
==21108==    by 0x4483E7: main (main.c:7)

I suspect this should be:

      matdat.datmat[i][4][0] = vec[0];
      matdat.datmat[i][4][1] = vec[1];
      matdat.datmat[i][4][2] = vec[2];

The same pattern of invalid access appears for the "Body Force Pressure Gradient":
https://github.com/PHASTA/phasta/blame/06d45c3e23cf9519431c2d5c7bbad33ddf2fbe7a/phSolver/common/input_fform.cc#L492-L494

Am I missing something here?

[hapfang]

I agree with Cameron, and if only two phases are simulated by using level set field, it then should be (index i should be 0 and 1, any explanation about why i=2 is appreciated.)

    for(i=0; i< levlset.iLSet ; i++){
      matdat.datmat[i][4][0] = vec[0];
      matdat.datmat[i][4][1] = vec[1];
      matdat.datmat[i][4][2] = vec[2];
    }

What concerns me is really how "Body Force" matrix datmat(1:2,5,1:3) is used in the source code. Take incompressible/e3res.f as an example:

phasta/phSolver/incompressible/e3res.f

Lines 378 to 380 in 06d45c3

	bfx = datmat(1,5,1) ! Boussinesq, g*alfap
	bfy = datmat(2,5,1)
	bfz = datmat(3,5,1)

Why is there a switch between 1st dimension and 3rd dimension of datmat? I would expect it to be something like

         bfx      = datmat(1,5,1) ! Boussinesq, g*alfap
         bfy      = datmat(1,5,2)
         bfz      = datmat(1,5,3)

Any clarification is appreciated.