Extractors for local_rhs in step-22

[Marco Feder]

Dear all,

I was going through step-22 and I saw that since you're using primitive elements, then you compute the local_rhs contribution using fe.system_to_component_index(i).first

However, as written in the program, we could as well multiply the dim+1 tensor having the values of the i-th shape function with the whole rhs vector (f,0). 

So, I have two questions: 

1) Are you using fe.system_to_component_index(i).first just to make the program run faster? Why not using extractors?

2) Taking the scalar product between a rank-1 Tensor and a Vector with the same size is not "conceptually" allowed. Did you have in mind something like the following?

//inside assemble_system()

std::vector<Tensor<1, dim>> phi_u(dofs_per_cell);

for (unsigned int q = 0; q < n_q_points; ++q)

{

for (unsigned int k = 0; k < dofs_per_cell; ++k)

{

//store symgrad_phi_u, div_phi_u, etc. using extractors

phi_u[k] = fe_values[velocities].value(k, q);

}

// compute local_matrix and local_preconditioner

double sum = 0.0; //store scalar product between first dim components of i-th shape function (corresponding to velocity) and first dim components of the rhs

for (unsigned int s = 0; s < dim; ++s)

{

sum += phi_u[i][s] * rhs_values[q][s];

}

local_rhs(i) +=

(sum + rhs_values[q][dim] * phi_p[i]) * fe_values.JxW(q);

}

}

or is there a better way to achieve this?

Best,
Marco

[Wolfgang Bangerth]

Marco,

> I was going through step-22 and I saw that since you're using primitive
> elements, then you compute the local_rhs contribution using

> /fe.system_to_component_index(i).first/

>
> However, as written in the program, we could as well multiply the dim+1 tensor
> having the values of the i-th shape function with the whole rhs vector (f,0).
>
> So, I have two questions:

> 1) Are you using /fe.system_to_component_index(i).first /just to make the

> program run faster? Why not using extractors?

I think it is an oversight. I would accept a patch to correct that :-)

> 2) Taking the scalar product between a rank-1 Tensor and a Vector with the
> same size is not "conceptually" allowed. Did you have in mind something like
> the following?
>
> //inside assemble_system()
> std::vector<Tensor<1, dim>> phi_u(dofs_per_cell);
> for (unsigned int q = 0; q < n_q_points; ++q)
> {
> for (unsigned int k = 0; k < dofs_per_cell; ++k)
> {
> //store symgrad_phi_u, div_phi_u, etc. using extractors
> phi_u[k] = fe_values[velocities].value(k, q);
> }
> // compute local_matrix and local_preconditioner
>
> double sum = 0.0; //store scalar product between first dim components of i-th
> shape function (corresponding to velocity) and first dim components of the rhs
> for (unsigned int s = 0; s < dim; ++s)
> {
> sum += phi_u[i][s] * rhs_values[q][s];
> }
> local_rhs(i) +=
> (sum + rhs_values[q][dim] * phi_p[i]) * fe_values.JxW(q);
> }
> }

That would work. I think a better solution would be to say outright that only
the first equation can have force terms, and then make the RightHandSide class
derived from TensorFunction<1,dim> and let it return its values as
Tensor<1,dim> objects for which you can take the dot product
fe_values[u].shape_value(q) * rhs_values[q]

Want to give this a try and write a patch? We'd be happy to walk you through
what is necessary for such a patch!

Best
Wolfgang

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Marco Feder]

Thanks Wolfgang,

indeed, as we have incompressibility, the last component of the rhs will always be 0. I wrote a possible patch, I'll open the pr tomorrow so that I can take another look at it!

Best,

Marco

[Wolfgang Bangerth]

On 11/19/21 6:26 PM, Marco Feder wrote:
>
> indeed, as we have incompressibility, the last component of the rhs will
> always be 0. I wrote a possible patch, I'll open the pr tomorrow so that I can
> take another look at it!

Awesome, thanks!
W.