Accessing vectorized arrays

[Patrick Greene]

Hello,

I have a basic question about accessing vectorized arrays, such as in the coefficient function in step 37. If I wanted to look at the z component of the points in the array and return a coefficient value of 1 when z>0 and a value of 2 when z<=0, for example, how would I do this? I have tried accessing them via p[2][0] and p[2][1], but this results in compiler errors saying "subscripted value is neither array nor pointer." Any advice would be greatly appreciated.

Thank you,

Patrick

[Martin Kronbichler]

Dear Patrick,

It is correct to access p[2][0] and p[2][1] if you want to access the z coordinates of the type Point<dim,VectorizedArray<double>>. But you get an error because step-37 channels through both Point<dim,double> (or just Point<dim> in most deal.II programs) and Point<dim,VectorizedArray<double> > through the same value method in the coefficient function, so then p[2][0] does not make sense for Point<dim,double>. To realize an access pattern as you describe it, you need to write two different functions for the two types separately.

Please let me know if you have any problems with that.

Best,

Martin 

[Patrick Greene]

Hi Martin,

Thank you for your reply. I made a template specialization for type double which seems to do the job. However, I am still having some trouble with outputting a vectorized array in the right form. I would like to write something like the following, which just checks the z-coordinate of the two points in the vectorized array and outputs a coefficient value for each of the points:

number coeff; //where "number" is VectorizedArray<double>

for (i=0; i<2; ++i) {

if (p[2][i] > 0)

coeff[i]=1.0

else coeff[i]=2.0

}

return coeff;

I get an error about index 1 not being in [0,1[, which I assume is due to not initializing the array correctly or not adding elements to it correctly, although I am not sure what the proper way is. Thanks again for any help.

[Martin Kronbichler]

Dear Patrick,

> number coeff; //where "number" is VectorizedArray<double>
> for (i=0; i<2; ++i) {
> if (p[2][i] > 0)
> coeff[i]=1.0
> else coeff[i]=2.0
> }
> return coeff;

It is generally unsafe to run loops from 0 to 2 like you write here. A
better form would be to run it up to the static member variable
VectorizedArray<number>::n_array_elements as the number depends on the
computer system. Your error is most likely because n_array_elements in
your vectorized array is 1 but you run the loop as if it were 2. Can you
check what VectorizedArray<number>::n_array_elements is in your case?
Cases currently implemented in deal.II have this value set to 1 except
for double & float number types where it is 2 & 4 for systems with SSE2
and 4 & 8 for systems with AVX.

> I get an error about index 1 not being in [0,1[, which I assume is due
> to not initializing the array correctly or not adding elements to it
> correctly, although I am not sure what the proper way is.

It shouldn't be because of improper initialization as VectorizedArray is
a POD data type where allocation = initialization (just as you have for
double, int, float and so on).

Best,
Martin

[Patrick Greene]

Hi Martin,

You were right, this was the problem. I assumed that vectorized arrays always had two entries, but the number of elements turned out to be one in my case. Thanks again for your help.

 -Patrick