Store data from a Tensor<1, dim> to a Vector<double>

[Judy Lee]

Hello everyone!

Can I create "Vector<double> x(n)", to store data values from a Tensor<1, dim> or Point<dim>, for physical coordinate after mesh/grid?

I can create "Point<dim> x", to store data values from a Tensor<1, dim>, it works thru like:

Point<dim> x;

for (const unsigned int i: fe_values.dof_indices())) {

x = cell->vertex(i);

}

std::cout << x << std::endl;

But how to use "Point<dim> x" for computing with other variables together? Like, I cannot implement "x * fe_values.shape_value(i, q_index)". However, it does not work thru if I tried like this: 

std::vector<Point<1>> x;

for (const unsigned int i: fe_values.dof_indices())) {

x[i] = cell->vertex(i);

}

std::cout << x[i] << std::endl;

It gives: 

make[3]: *** [CMakeFiles/run.dir/build.make:77: CMakeFiles/run] Segmentation fault

make[2]: *** [CMakeFiles/MakeFile2:237: CMakeFiles/run.dir/all] Error 2

make[1]: *** [CMakeFiles/MakeFile2:244: CMakeFiles/run.dir/rule] Error 2

make: *** [Makefile:202: run] Error 2

How to store data values from Tensor<1, dim>, into a created "Vector<double> x(n)", for computing with other variables together?

Thank you so much!

Judy

[Wolfgang Bangerth]

On 12/7/21 03:30, Judy Lee wrote:
>
> But how to use "Point<dim> x" for computing with other variables
> together? Like, I cannot implement "x * fe_values.shape_value(i,
> q_index)". However, it does not work thru if I tried like this:
>
> std::vector<Point<1>> x;
> for (const unsigned int i: fe_values.dof_indices())) {
> x[i] = cell->vertex(i);
> }
> std::cout << x[i] << std::endl;
>
> It gives:
> make[3]: *** [CMakeFiles/run.dir/build.make:77: CMakeFiles/run]
> Segmentation fault
> make[2]: *** [CMakeFiles/MakeFile2:237: CMakeFiles/run.dir/all] Error 2
> make[1]: *** [CMakeFiles/MakeFile2:244: CMakeFiles/run.dir/rule] Error 2
> make: *** [Makefile:202: run] Error 2

Judy -- can you explain *what* it is that you want to do? It is often
easier to answer *what* questions, but you ask a *how* question.

Separately, the error you show is a segmentation fault. I don't think it
has anything to do with the line

x[i] = cell->vertex(i)

but instead with the fact that your vector 'x' has size zero. You need
to set it to the correct size before you access elements 'x[i]'.

Best
W.

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

[Matthias Maier]

Dear Judy,

On Tue, Dec 7, 2021, at 11:45 CST, Judy Lee <reader....@gmail.com> wrote:

> cell_x(i) += fe_values.shape_value(i, q_index) * x(i);

Here, you are multiplying a scalar (fe_values.shape_value(i, q_index))
and a rank-1 tensor (x(i)), which results in a rank-1 tensor.

The error message says that you cannot store this rank-1 tensor in
cell_x(i), which is a scalar (vector<double>).

I don't know whether this is what you want to do, but in 1D you can
simply take the first value of your rank-1 tensor to transform the
expression into a scalar:

cell_x(i) += fe_values.shape_value(i, q_index) * x(i)[0];

However, this is not a good "dimension independent" approach, meaning
the moment you try to solve in 2D or 3D you will have to modify your
code.

Best,
Matthias

[Judy Lee]

Dear Matthias,

Thanks so much for your reply.

Your comment and Dr. Bangerth's comment both exactly pointed out the problems in my settings on the scalar "Vector<double> cell_x(dofs_per_cell)", and the rank-1 tensor "std::vector<Point<dim>> x(dofs_per_cell)", which multiply the scalar "fe_values.shape_value(i, q_index)" in the expression.

Yes, vertices of each cell form a rank-1 tensor. I just put in the first value of the rank-1 tensor, to transform the expression into the scalar "cell_x", that perfectly works thru. This is amazing!

Thank you all so much! I will keep up my work.

Best,
Judy

[Judy Lee]

Good afternoon, everyone!

I have a couple of follow-up questions.
(1) Is there an internal command/function that outputs higher-order-elemental (Q2, Q3 elements) nodal physical coordinate?
Just like, "cell->vertex" outputs grid points that works for Q1 element.
If I want to count in a linear function f(x) = x on Laplace's right-hand-side, should I just manually make an expression, to interpolate mid-nodes, as well as other fractional-nodes @ 1/3 and 2/3 of each cell "he", to construct “x”?

(2) In the nonlinear problem, are nodal physical coordinates updated thru finite element program, like, including "phi_i * x_i" or “grad_phi_i * x_i” for any i of each cell?

Thank you very much!
Best,
Judy

On Tuesday, December 7, 2021 at 1:44:38 PM UTC-5 Matthias Maier wrote:

[Wolfgang Bangerth]

Judy:

> (1) Is there an internal command/function that outputs
> higher-order-elemental (Q2, Q3 elements) nodal physical coordinate?
> Just like, "cell->vertex" outputs grid points that works for Q1 element.
> If I want to count in a linear function f(x) = x on Laplace's
> right-hand-side, should I just manually make an expression, to
> interpolate mid-nodes, as well as other fractional-nodes @ 1/3 and 2/3
> of each cell "he", to construct “x”?

You can query the location of all degrees of freedom via the
DoFTools::map_dofs_to_support_points() function. But I don't understand
the connection to your right hand side: you generally only ever need to
evaluate the right hand side function at *quadrature points*, not at
*node points*. You get the physical locations of quadrature points via
fe_values.quadrature_point(q_point)
as shown in several of the tutorial programs.

> (2) In the nonlinear problem, are nodal physical coordinates updated
> thru finite element program, like, including "phi_i * x_i" or
> “grad_phi_i * x_i” for any i of each cell?

The library doesn't do anything unless you instruct it to, because it
has no idea what kind of equation you are solving. So, unless you
explicitly say "update the node location with this kind of information",
nothing will happen.

[Judy Lee]

Dr. Bangerth,

Thank you so much!

I will keep up my work.

Best,
Judy

[Wolfgang Bangerth]

On 12/14/21 08:32, Judy Lee wrote:
>
> I tried *DoFTools::map_dofs_to_support_points<dim>(mapping, dof_handler,
> support_points)*, to obtain coordinates on *local_dof_indices*, not
> running thru yet.
> My .cc file is attached with reference to step-3 in 2D, in which I added
> Lines 137-138 with reference to Lines 126-127, as well as by looking at
> step-2 introduction.
> Error shows, " ’mapping’ was not declared in the scope". I added Line 62
> that does not figure out the error, with adding Line 40.

If you do this, you get the following error message:

/home/fac/g/bangerth/p/deal.II/1/install/examples/step-3/step-3.cc:67:30:
error: no matching function for call to ‘dealii::MappingQ<2, 2>::MappingQ()’
, dof_handler(triangulation)

It tells you that you need to explicitly initialize the 'mapping'
variable in the constructor of the class, in the same way as you do for
the finite element.

[Judy Lee]

Dr. Bangerth,

Thank you so much!

Points are printed, I will keep up my work.

Best,
Judy