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

Il giorno venerdì 19 novembre 2021 alle 23:15:50 UTC+1 Wolfgang Bangerth ha 
scritto:

>
> 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: [email protected]
> www: http://www.math.colostate.edu/~bangerth/
>
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to [email protected].
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/d83ed03f-5efd-4d94-a97b-b250daa51954n%40googlegroups.com.

Reply via email to