[deal.II] Re: Vector-valued gradient of solution vector

[Daniel Arndt]

Joel,

Yes, this looks about right. I would really replace "std::vector double 
x,y,z" by "Vector x,y,z".
In the end, entries in these vectors only make sense with the corresponding 
DoFHandler.
This would also allow you to use these Vectors for postprocessing in 
deal.II.

Best,
Daniel

-- 
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 dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Vector-valued gradient of solution vector

[Joel Davidsson]

Dear Daniel,

Ok, I think I got it to work now, can you verify that it is correct or if I 
have missed something crucial.

  std::vector   x,y,z;
  x.resize (dof_handler_scalar.n_dofs());
  y.resize (dof_handler_scalar.n_dofs());
  z.resize (dof_handler_scalar.n_dofs());

 Quadrature quadrature_formula(fe_scalar.get_unit_support_points());

 FEValues fe_values (fe, quadrature_formula,
  update_values | update_quadrature_points);

 FEValues fe_values_scalar (fe_scalar, quadrature_formula,
  update_values | update_quadrature_points);

 std::vector> local_velocity_values 
(quadrature_formula.size());

 for (unsigned int i=0;i local_dof_indices (dofs_per_cell);

 typename DoFHandler::active_cell_iterator
   cell_scalar = dof_handler_scalar.begin_active(),
   endc_scalar = dof_handler_scalar.end();

 typename DoFHandler::active_cell_iterator
 cell = dof_handler.begin_active(),
 endc = dof_handler.end();

 for (; cell!=endc; ++cell, ++cell_scalar)
   {
 fe_values.reinit (cell);
 fe_values_scalar.reinit (cell_scalar);
 fe_values.get_function_values (solution,local_velocity_values);

 cell_scalar->get_dof_indices (local_dof_indices);

   for (unsigned int i=0; i
> Joel,
>
> I did try it, but I cant just change the local_velocity_values from a 
>> std::vector> to std::vector. Then the 
>> fe_values[velocities].get_function_values function doesnt compile. 
>>
> That's true, but you should use std::vector > and not 
> std::vector, see also [1].
>  
>
>> (correct me if wrong) This method of using norm_sqr() would only produce 
>> the size, not the individual x,y and z components. I need the x,y and 
>> z components separately so I can get the size but also the difference 
>> between two solutions and that error of these two solutions. Hence my 
>> previous message. I want to be able to manipulate the vector components. 
>> How would I do this?
>>
> Yes, but since you get all the values at the local dofs, you can play 
> around with them as you like and store the result in velocity. If this is 
> not sufficient, you can also make velocity a std::vector and use something 
> like
> velocity[component][local_dof_indices[i]] = 
> local_velocity_values[i][component]
>
> Best,
> Daniel
>
>>
> [1] 
> https://www.dealii.org/8.4.1/doxygen/deal.II/classFEValuesBase.html#a396f1430dd3f5716a9fe6ef2762edb5d
>

-- 
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 dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Vector-valued gradient of solution vector

[Daniel Arndt]

Joel,

I did try it, but I cant just change the local_velocity_values from a 
> std::vector> to std::vector. Then the 
> fe_values[velocities].get_function_values function doesnt compile. 
>
That's true, but you should use std::vector > and not 
std::vector, see also [1].
 

> (correct me if wrong) This method of using norm_sqr() would only produce 
> the size, not the individual x,y and z components. I need the x,y and 
> z components separately so I can get the size but also the difference 
> between two solutions and that error of these two solutions. Hence my 
> previous message. I want to be able to manipulate the vector components. 
> How would I do this?
>
Yes, but since you get all the values at the local dofs, you can play 
around with them as you like and store the result in velocity. If this is 
not sufficient, you can also make velocity a std::vector and use something 
like
velocity[component][local_dof_indices[i]] = 
local_velocity_values[i][component]

Best,
Daniel

>
[1] 
https://www.dealii.org/8.4.1/doxygen/deal.II/classFEValuesBase.html#a396f1430dd3f5716a9fe6ef2762edb5d

-- 
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 dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Vector-valued gradient of solution vector

[Joel Davidsson]

Dear Daniel,

I did try it, but I cant just change the local_velocity_values from a 
std::vector> to std::vector. Then the 
fe_values[velocities].get_function_values function doesnt compile. So I 
commented that part out, the results you see is from the code as it is, I 
used the for loop at line 354 to get the size (which I know is wrong, but 
at least it give something).  Even if I had found the correct way to do 
this, it would not have solved my original problem.

(correct me if wrong) This method of using norm_sqr() would only produce 
the size, not the individual x,y and z components. I need the x,y and 
z components separately so I can get the size but also the difference 
between two solutions and that error of these two solutions. Hence my 
previous message. I want to be able to manipulate the vector components. 
How would I do this?

Best,
Joel


On Tuesday, September 27, 2016 at 7:38:48 PM UTC+2, Daniel Arndt wrote:
>
> Joel,
>
> Its been a week since my latest post, I was just wondering if you could 
>> help me.
>>
> Did you try what I suggested? What were the problems?
>
> Looking at the latest files you sent, local_velocity_values has still not 
> the type std::vector >
> and you are still using
>
> for (unsigned int q_index=0; q_index   for (unsigned int i=0; i velocity[local_dof_indices[i]] = local_velocity_values[q_index];
>
> instead of
>
> for (unsigned int i=0; i   velocity[local_dof[indices[i]] = local_velocity_values(i).norm_sqr();
>
> This way you don't rely on any particular order of the dofs.
>
> Best,
> Daniel
>

-- 
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 dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Vector-valued gradient of solution vector

[Daniel Arndt]

Joel,

Its been a week since my latest post, I was just wondering if you could 
> help me.
>
Did you try what I suggested? What were the problems?

Looking at the latest files you sent, local_velocity_values has still not 
the type std::vector >
and you are still using

for (unsigned int q_index=0; q_indexhttp://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 dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Vector-valued gradient of solution vector

[Joel Davidsson]

Dear all,

Its been a week since my latest post, I was just wondering if you could 
help me.

Best,
Joel

On Tuesday, September 20, 2016 at 9:12:23 AM UTC+2, Joel Davidsson wrote:
>
> Dear Daniel,
>
> To clarify what I am searching for, I would like to get the vector-valued 
> solution into a format that matches a scalar. Either as a 
> std::vector> or a std::vector>, the format doesn't 
> matter as long as it matches a scalar dofs. So I can manipulate the vector 
> components and for example get the size.
>
> I made a test program that I have attached, that takes the gradient of a 
> gaussian function. At line 354, I do the convergence from a vector valued 
> solution to get the size of the vector (as I posted before). This was wrong 
> and you guys said that I shouldn't assume any particular order of the dofs. 
> Then I tried to this conversion the correct way, see line 365, but I do not 
> know how to do this conversion the correct way.
>
> The code I have sent you, produces the correct result for fe1 and fe2, 
> i.e. when I set fe_scalar (2), fe(FE_Q(2),dim) at line 100. But failed 
> for fe3. I have attached the result I got for fe2 and fe3 as well as the 
> plot script I have used to generate these plots.
>
> I would be really grateful if you could help me figure out how to solve 
> this.
>
> Best,
> Joel
>
> On Monday, September 19, 2016 at 11:56:51 AM UTC+2, Daniel Arndt wrote:
>>
>> Joel,
>>
>>> [...]
>>>  for (unsigned int q_index=0; q_index>>for (unsigned int i=0; i>>  {
>>> velocity[local_dof_indices[i]] = local_velocity_values[q_index];
>>>  }
>>>
>>> Here you assign the same Tensor<1,dim> at all DoFs multiple times and 
>> you end with velocity[*]=local_velocity_values[n_q_points]-1.
>> I don't think that this is what you want to do, is it? If I understand 
>> you correctly, you try to assign the magnitude of the velocity in each 
>> degree of freedom
>> to a different scalar Vector that represents a scalar FiniteElement 
>> Vector.
>> In this case, you need to use a Quadrature object that uses the 
>> unit_support_points of the scalar FE [1] and your assignment would read
>>
>> for (unsigned int i=0; i>   velocity[local_dof[indices[i]] = local_velocity_values(i).norm_sqr();
>>
>> where velocity is a Vector that is related to fe_scalar and 
>> local_velocity_values has the type  std::vector >.
>>
>> Best,
>> Daniel
>>
>> [1] 
>> https://github.com/dealii/dealii/wiki/Frequently-Asked-Questions#how-to-get-the-mapped-position-of-support-points-of-my-element
>>
>

-- 
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 dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Vector-valued gradient of solution vector

[Joel Davidsson]

Dear Daniel,

To clarify what I am searching for, I would like to get the vector-valued 
solution into a format that matches a scalar. Either as a 
std::vector> or a std::vector>, the format doesn't 
matter as long as it matches a scalar dofs. So I can manipulate the vector 
components and for example get the size.

I made a test program that I have attached, that takes the gradient of a 
gaussian function. At line 354, I do the convergence from a vector valued 
solution to get the size of the vector (as I posted before). This was wrong 
and you guys said that I shouldn't assume any particular order of the dofs. 
Then I tried to this conversion the correct way, see line 365, but I do not 
know how to do this conversion the correct way.

The code I have sent you, produces the correct result for fe1 and fe2, i.e. 
when I set fe_scalar (2), fe(FE_Q(2),dim) at line 100. But failed for 
fe3. I have attached the result I got for fe2 and fe3 as well as the plot 
script I have used to generate these plots.

I would be really grateful if you could help me figure out how to solve 
this.

Best,
Joel

On Monday, September 19, 2016 at 11:56:51 AM UTC+2, Daniel Arndt wrote:
>
> Joel,
>
>> [...]
>>  for (unsigned int q_index=0; q_index>for (unsigned int i=0; i>  {
>> velocity[local_dof_indices[i]] = local_velocity_values[q_index];
>>  }
>>
>> Here you assign the same Tensor<1,dim> at all DoFs multiple times and you 
> end with velocity[*]=local_velocity_values[n_q_points]-1.
> I don't think that this is what you want to do, is it? If I understand you 
> correctly, you try to assign the magnitude of the velocity in each degree 
> of freedom
> to a different scalar Vector that represents a scalar FiniteElement Vector.
> In this case, you need to use a Quadrature object that uses the 
> unit_support_points of the scalar FE [1] and your assignment would read
>
> for (unsigned int i=0; i   velocity[local_dof[indices[i]] = local_velocity_values(i).norm_sqr();
>
> where velocity is a Vector that is related to fe_scalar and 
> local_velocity_values has the type  std::vector >.
>
> Best,
> Daniel
>
> [1] 
> https://github.com/dealii/dealii/wiki/Frequently-Asked-Questions#how-to-get-the-mapped-position-of-support-points-of-my-element
>

-- 
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 dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
/* -
 *
 * Copyright (C) 1999 - 2013 by the deal.II authors
 *
 * This file is part of the deal.II library.
 *
 * The deal.II library is free software; you can use it, redistribute
 * it, and/or modify it under the terms of the GNU Lesser General
 * Public License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 * The full text of the license can be found in the file LICENSE at
 * the top level of the deal.II distribution.
 *
 * -

 *
 * Author: Wolfgang Bangerth, University of Heidelberg, 1999
 */

#include 

#include 
#include 
#include 

#include 
#include 
#include 
#include 

#include 
#include 
#include 
#include 

#include 
#include 

#include 
#include 

#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 

#include 
#include 
#include 
#include 

#include 

using namespace dealii;

const double pi=3.14159265359;
unsigned int center_id=0;

template 
class Step4
{
public:
  Step4 ();
  void run ();

private:
  void make_grid ();
  void setup_system();
  void assemble_system ();
  void solve ();
  void output_results () const;

  Triangulation   triangulation;
  FE_Qfe_scalar;
  FESystem		   fe;

  DoFHandler  dof_handler_scalar;
  DoFHandler  dof_handler;

  SparsityPattern  sparsity_pattern;
  SparseMatrix system_matrix;

  ConstraintMatrix constraints;

  Vector   solution;
  Vector   system_rhs;
  Vector   test;
  Vector   size;

  std::vector> nodes;
};

template 
Step4::Step4 ()
  :
  fe_scalar (3),
  fe (FE_Q(3),dim),
  dof_handler (triangulation),
  dof_handler_scalar (triangulation)
{}

template 
void Step4::make_grid ()
{
  GridGenerator::hyper_cube (triangulation, -2.5, 2.5);
  triangulation.refine_global (2);

  for (unsigned int step=0; step<2; ++step)
  {
	  typename Triangulation::active_cell_iterator
	  cell = triangulation.begin_active(),
	  endc = triangulation.end();
	  for (; cell!=endc; ++cell)
	  {
		  for (unsigned int v=0;v < GeometryInfo::vertices_per_cell;++v)
		  {
			  Point center;
			  for (unsigned int i=0; ivertex(v));
			 

[deal.II] Re: Vector-valued gradient of solution vector

[Daniel Arndt]

Joel,

> [...]
>  for (unsigned int q_index=0; q_indexfor (unsigned int i=0; i  {
> velocity[local_dof_indices[i]] = local_velocity_values[q_index];
>  }
>
> Here you assign the same Tensor<1,dim> at all DoFs multiple times and you 
end with velocity[*]=local_velocity_values[n_q_points]-1.
I don't think that this is what you want to do, is it? If I understand you 
correctly, you try to assign the magnitude of the velocity in each degree 
of freedom
to a different scalar Vector that represents a scalar FiniteElement Vector.
In this case, you need to use a Quadrature object that uses the 
unit_support_points of the scalar FE [1] and your assignment would read

for (unsigned int i=0; i >.

Best,
Daniel

[1] 
https://github.com/dealii/dealii/wiki/Frequently-Asked-Questions#how-to-get-the-mapped-position-of-support-points-of-my-element

-- 
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 dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Vector-valued gradient of solution vector

[Joel Davidsson]

Dear J-P,

Ok, I tried to to use the FEValueExtractors but I am not sure if my code is 
working correctly. 

   std::vector>   velocity;
   velocity.resize (dof_handler_scalar.n_dofs());

 QGauss  quadrature_formula(dof_handler.get_fe().degree+1);

 FEValues fe_values (fe, quadrature_formula,
  update_values   | update_gradients |
  update_quadrature_points | update_JxW_values);

 FEValues fe_values_scalar (fe_scalar, quadrature_formula,
  update_values   | update_gradients |
  update_quadrature_points | update_JxW_values);

 const unsigned int   dofs_per_cell = fe_scalar.dofs_per_cell;
 const unsigned int   n_q_points= quadrature_formula.size();

 std::vector local_dof_indices (dofs_per_cell);
  std::vector > local_velocity_values (n_q_points);

 const FEValuesExtractors::Vector velocities (0);

 typename DoFHandler::active_cell_iterator
   cell_scalar = dof_handler_scalar.begin_active(),
   endc_scalar = dof_handler_scalar.end();

 typename DoFHandler::active_cell_iterator
 cell = dof_handler.begin_active(),
 endc = dof_handler.end();

 for (; cell!=endc; ++cell, ++cell_scalar)
   {
 fe_values.reinit (cell);
 fe_values_scalar.reinit (cell_scalar);
 fe_values[velocities].get_function_values 
(solution,local_velocity_values);

 cell_scalar->get_dof_indices (local_dof_indices);

 for (unsigned int q_index=0; q_index
> Dear Joel,
>
> In general, no you cannot. Here you assume a specific ordering of the DoFs 
> in the global solution vector, and this changes for each finite element. 
> For example, using FE_Q's might result the ordering [x0 y0 z0 x1 y1 z1 ...] 
> and FE_DGQ's might lead to [x0 x1 x2 ... y0 y1 y2 ... z0 z1 z2...]. To do 
> this safely you should definitely do this using the built in tools, such as 
> FEValueExtractors 
> 
>  
> or  fe.system_to_component_index 
> .
>  
> See the tutorials for examples on how to use them.
>
> Regards,
> J-P
>
> On Tuesday, September 13, 2016 at 2:05:56 PM UTC+2, Joel Davidsson wrote:
>>
>> Dear all,
>>
>> If I have a Vector-valued solution in 3D and I would like to get the size 
>> of each vector in another smaller Vector, can I do it like this:
>>
>> for (unsigned int i=0;i> {
>>  unsigned int x=i;
>>  unsigned int y=i+1;
>>  unsigned int z=i+2;
>>  unsigned int j = i / 3.0;
>>
>>  size(j)=pow(solution(x),2.0)+pow(solution(y),2.0)+pow(solution(z),2.0);
>> }
>>
>> Where solution is the vector valued Vector that is 3 times bigger that 
>> the scalar Vector called size. 
>>
>> Best,
>> Joel
>>
>> On Tuesday, August 30, 2016 at 6:48:46 PM UTC+2, Jean-Paul Pelteret wrote:
>>>
>>> Dear Joel,
>>>
>>> No, this functionality doesn't exist. However, you can create a 
>>> quadrature rule with the quadrature points located at the support points 
>>> (see the FAQ 
>>> ),
>>>  
>>> and could then extract the solution based on your custom quadrature rule.
>>>
>>> Regards,
>>> J-P
>>>
>>> On Tuesday, August 30, 2016 at 5:32:13 PM UTC+2, Joel Davidsson wrote:

 Dear Daniel,

 I was just wondering if their exist a similar function 
 to DoFTools::map_dofs_to_support_points, but that maps a vector-valued 
 solution Vector into to a std::vector> instead?

 Best,
 Joel

 On Monday, August 22, 2016 at 9:03:22 PM UTC+2, Daniel Arndt wrote:
>
> Joel,
>
> You should not rely on any particular sorting of the dofs in the 
> solution vector. Instead you can ask the FiniteElement on each cell for 
> the 
> component of a local dof by
> const unsigned int component = fe.system_to_component_index(i).first; 
>
> Apart from that the DataOut object knows about components and you can 
> do some postprocessing on the solution vector using DataPostprocessor [1]
>
> Best,
> Daniel
>
>  [1] 
> https://dealii.org/8.4.1/doxygen/deal.II/classDataPostprocessor.html
>
> Am Montag, 22. August 2016 17:05:08 UTC+2 schrieb Joel Davidsson:
>>
>> Dear Daniel,
>>
>> Thank you for a very good answer, adding the fe_values_scalar and 
>> cell_scalar fixed the problem.
>>
>> I have a follow-up question about the solution I get out, how is the 
>> data organized in the solution vector? Say for example I want to loop 
>> over 
>> all the x components, how would I do that?
>>
>> Best,
>> Joel
>>
>> On Friday, August 12, 2016 at 5:35:07 PM UTC+2, Daniel Arndt wrote:
>>>
>>> Joel,
>>>
>>> Yes the matrix you are assembling is a vector-valued mass matrix now.
>>> For me your

[deal.II] Re: Vector-valued gradient of solution vector

[Jean-Paul Pelteret]

Dear Joel,

In general, no you cannot. Here you assume a specific ordering of the DoFs 
in the global solution vector, and this changes for each finite element. 
For example, using FE_Q's might result the ordering [x0 y0 z0 x1 y1 z1 ...] 
and FE_DGQ's might lead to [x0 x1 x2 ... y0 y1 y2 ... z0 z1 z2...]. To do 
this safely you should definitely do this using the built in tools, such as 
FEValueExtractors 

 
or  fe.system_to_component_index 
.
 
See the tutorials for examples on how to use them.

Regards,
J-P

On Tuesday, September 13, 2016 at 2:05:56 PM UTC+2, Joel Davidsson wrote:
>
> Dear all,
>
> If I have a Vector-valued solution in 3D and I would like to get the size 
> of each vector in another smaller Vector, can I do it like this:
>
> for (unsigned int i=0;i {
>  unsigned int x=i;
>  unsigned int y=i+1;
>  unsigned int z=i+2;
>  unsigned int j = i / 3.0;
>
>  size(j)=pow(solution(x),2.0)+pow(solution(y),2.0)+pow(solution(z),2.0);
> }
>
> Where solution is the vector valued Vector that is 3 times bigger that the 
> scalar Vector called size. 
>
> Best,
> Joel
>
> On Tuesday, August 30, 2016 at 6:48:46 PM UTC+2, Jean-Paul Pelteret wrote:
>>
>> Dear Joel,
>>
>> No, this functionality doesn't exist. However, you can create a 
>> quadrature rule with the quadrature points located at the support points 
>> (see the FAQ 
>> ),
>>  
>> and could then extract the solution based on your custom quadrature rule.
>>
>> Regards,
>> J-P
>>
>> On Tuesday, August 30, 2016 at 5:32:13 PM UTC+2, Joel Davidsson wrote:
>>>
>>> Dear Daniel,
>>>
>>> I was just wondering if their exist a similar function 
>>> to DoFTools::map_dofs_to_support_points, but that maps a vector-valued 
>>> solution Vector into to a std::vector> instead?
>>>
>>> Best,
>>> Joel
>>>
>>> On Monday, August 22, 2016 at 9:03:22 PM UTC+2, Daniel Arndt wrote:

 Joel,

 You should not rely on any particular sorting of the dofs in the 
 solution vector. Instead you can ask the FiniteElement on each cell for 
 the 
 component of a local dof by
 const unsigned int component = fe.system_to_component_index(i).first; 

 Apart from that the DataOut object knows about components and you can 
 do some postprocessing on the solution vector using DataPostprocessor [1]

 Best,
 Daniel

  [1] 
 https://dealii.org/8.4.1/doxygen/deal.II/classDataPostprocessor.html

 Am Montag, 22. August 2016 17:05:08 UTC+2 schrieb Joel Davidsson:
>
> Dear Daniel,
>
> Thank you for a very good answer, adding the fe_values_scalar and 
> cell_scalar fixed the problem.
>
> I have a follow-up question about the solution I get out, how is the 
> data organized in the solution vector? Say for example I want to loop 
> over 
> all the x components, how would I do that?
>
> Best,
> Joel
>
> On Friday, August 12, 2016 at 5:35:07 PM UTC+2, Daniel Arndt wrote:
>>
>> Joel,
>>
>> Yes the matrix you are assembling is a vector-valued mass matrix now.
>> For me your code is failing with
>>
>> void dealii::FEValuesBase> spacedim>::get_function_gradients(const InputVector&, 
>> std::vector> InputVector::value_type> 
>> >&) const [with InputVector = dealii::Vector; int dim = 3; int 
>> spacedim = 3; typename InputVector::value_type = double]
>> The violated condition was: 
>> (fe->n_components()) == (1)
>> The name and call sequence of the exception was:
>> dealii::ExcDimensionMismatch((fe->n_components()),(1))
>> Additional Information: 
>> Dimension 3 not equal to 1
>>
>> and this is not surprising. What you are doing is to extract the 
>> gradients of a vector-valued finite element solution. The object that 
>> should store this should therefore be a 
>> std::vector> as you want to store for each 
>> quadrature point and each component a Tensor<1,dim>.
>>
>> What you want to do is really that in "Do not work". As you have two 
>> DoFHandlers, you should also have two FEValues objects and two 
>> cell_iterator corresponding to the correct DoFHandler. Then you can 
>> extract 
>> the gradient of your scalar field and project it onto the ansatz space 
>> for 
>> the vector-valued field. This should look like:
>>
>> typename DoFHandler::active_cell_iterator cell_scalar = 
>> dof_handler_scalar.begin_active();
>> typename DoFHandler::active_cell_iterator cell = 
>> dof_handler.begin_active();
>> ...
>>
>>  for (; cell!=endc; ++cell, ++cell_vector)
>> {
>>   fe_values.reinit (ce

[deal.II] Re: Vector-valued gradient of solution vector

[Joel Davidsson]

Dear all,

If I have a Vector-valued solution in 3D and I would like to get the size 
of each vector in another smaller Vector, can I do it like this:

for (unsigned int i=0;i
> Dear Joel,
>
> No, this functionality doesn't exist. However, you can create a quadrature 
> rule with the quadrature points located at the support points (see the FAQ 
> ),
>  
> and could then extract the solution based on your custom quadrature rule.
>
> Regards,
> J-P
>
> On Tuesday, August 30, 2016 at 5:32:13 PM UTC+2, Joel Davidsson wrote:
>>
>> Dear Daniel,
>>
>> I was just wondering if their exist a similar function 
>> to DoFTools::map_dofs_to_support_points, but that maps a vector-valued 
>> solution Vector into to a std::vector> instead?
>>
>> Best,
>> Joel
>>
>> On Monday, August 22, 2016 at 9:03:22 PM UTC+2, Daniel Arndt wrote:
>>>
>>> Joel,
>>>
>>> You should not rely on any particular sorting of the dofs in the 
>>> solution vector. Instead you can ask the FiniteElement on each cell for the 
>>> component of a local dof by
>>> const unsigned int component = fe.system_to_component_index(i).first; 
>>>
>>> Apart from that the DataOut object knows about components and you can do 
>>> some postprocessing on the solution vector using DataPostprocessor [1]
>>>
>>> Best,
>>> Daniel
>>>
>>>  [1] 
>>> https://dealii.org/8.4.1/doxygen/deal.II/classDataPostprocessor.html
>>>
>>> Am Montag, 22. August 2016 17:05:08 UTC+2 schrieb Joel Davidsson:

 Dear Daniel,

 Thank you for a very good answer, adding the fe_values_scalar and 
 cell_scalar fixed the problem.

 I have a follow-up question about the solution I get out, how is the 
 data organized in the solution vector? Say for example I want to loop over 
 all the x components, how would I do that?

 Best,
 Joel

 On Friday, August 12, 2016 at 5:35:07 PM UTC+2, Daniel Arndt wrote:
>
> Joel,
>
> Yes the matrix you are assembling is a vector-valued mass matrix now.
> For me your code is failing with
>
> void dealii::FEValuesBase::get_function_gradients(const 
> InputVector&, std::vector InputVector::value_type> >&) const [with InputVector = 
> dealii::Vector; int dim = 3; int spacedim = 3; typename 
> InputVector::value_type = double]
> The violated condition was: 
> (fe->n_components()) == (1)
> The name and call sequence of the exception was:
> dealii::ExcDimensionMismatch((fe->n_components()),(1))
> Additional Information: 
> Dimension 3 not equal to 1
>
> and this is not surprising. What you are doing is to extract the 
> gradients of a vector-valued finite element solution. The object that 
> should store this should therefore be a 
> std::vector> as you want to store for each 
> quadrature point and each component a Tensor<1,dim>.
>
> What you want to do is really that in "Do not work". As you have two 
> DoFHandlers, you should also have two FEValues objects and two 
> cell_iterator corresponding to the correct DoFHandler. Then you can 
> extract 
> the gradient of your scalar field and project it onto the ansatz space 
> for 
> the vector-valued field. This should look like:
>
> typename DoFHandler::active_cell_iterator cell_scalar = 
> dof_handler_scalar.begin_active();
> typename DoFHandler::active_cell_iterator cell = 
> dof_handler.begin_active();
> ...
>
>  for (; cell!=endc; ++cell, ++cell_vector)
> {
>   fe_values.reinit (cell);
>   fe_values_scalar.reinit (cell_scalar);
>   fe_values_scalar.get_function_gradients(test,fe_function_grad);
>   ...
> }
>
> Best,
> Daniel
>


-- 
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 dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Vector-valued gradient of solution vector

[Jean-Paul Pelteret]

Dear Joel,

No, this functionality doesn't exist. However, you can create a quadrature 
rule with the quadrature points located at the support points (see the FAQ 
),
 
and could then extract the solution based on your custom quadrature rule.

Regards,
J-P

On Tuesday, August 30, 2016 at 5:32:13 PM UTC+2, Joel Davidsson wrote:
>
> Dear Daniel,
>
> I was just wondering if their exist a similar function 
> to DoFTools::map_dofs_to_support_points, but that maps a vector-valued 
> solution Vector into to a std::vector> instead?
>
> Best,
> Joel
>
> On Monday, August 22, 2016 at 9:03:22 PM UTC+2, Daniel Arndt wrote:
>>
>> Joel,
>>
>> You should not rely on any particular sorting of the dofs in the solution 
>> vector. Instead you can ask the FiniteElement on each cell for the 
>> component of a local dof by
>> const unsigned int component = fe.system_to_component_index(i).first; 
>>
>> Apart from that the DataOut object knows about components and you can do 
>> some postprocessing on the solution vector using DataPostprocessor [1]
>>
>> Best,
>> Daniel
>>
>>  [1] https://dealii.org/8.4.1/doxygen/deal.II/classDataPostprocessor.html
>>
>> Am Montag, 22. August 2016 17:05:08 UTC+2 schrieb Joel Davidsson:
>>>
>>> Dear Daniel,
>>>
>>> Thank you for a very good answer, adding the fe_values_scalar and 
>>> cell_scalar fixed the problem.
>>>
>>> I have a follow-up question about the solution I get out, how is the 
>>> data organized in the solution vector? Say for example I want to loop over 
>>> all the x components, how would I do that?
>>>
>>> Best,
>>> Joel
>>>
>>> On Friday, August 12, 2016 at 5:35:07 PM UTC+2, Daniel Arndt wrote:

 Joel,

 Yes the matrix you are assembling is a vector-valued mass matrix now.
 For me your code is failing with

 void dealii::FEValuesBase::get_function_gradients(const 
 InputVector&, std::vector>>> InputVector::value_type> >&) const [with InputVector = 
 dealii::Vector; int dim = 3; int spacedim = 3; typename 
 InputVector::value_type = double]
 The violated condition was: 
 (fe->n_components()) == (1)
 The name and call sequence of the exception was:
 dealii::ExcDimensionMismatch((fe->n_components()),(1))
 Additional Information: 
 Dimension 3 not equal to 1

 and this is not surprising. What you are doing is to extract the 
 gradients of a vector-valued finite element solution. The object that 
 should store this should therefore be a 
 std::vector> as you want to store for each 
 quadrature point and each component a Tensor<1,dim>.

 What you want to do is really that in "Do not work". As you have two 
 DoFHandlers, you should also have two FEValues objects and two 
 cell_iterator corresponding to the correct DoFHandler. Then you can 
 extract 
 the gradient of your scalar field and project it onto the ansatz space for 
 the vector-valued field. This should look like:

 typename DoFHandler::active_cell_iterator cell_scalar = 
 dof_handler_scalar.begin_active();
 typename DoFHandler::active_cell_iterator cell = 
 dof_handler.begin_active();
 ...

  for (; cell!=endc; ++cell, ++cell_vector)
 {
   fe_values.reinit (cell);
   fe_values_scalar.reinit (cell_scalar);
   fe_values_scalar.get_function_gradients(test,fe_function_grad);
   ...
 }

 Best,
 Daniel

>>>

-- 
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 dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Vector-valued gradient of solution vector

[Joel Davidsson]

Dear Daniel,

I was just wondering if their exist a similar function 
to DoFTools::map_dofs_to_support_points, but that maps a vector-valued 
solution Vector into to a std::vector> instead?

Best,
Joel

On Monday, August 22, 2016 at 9:03:22 PM UTC+2, Daniel Arndt wrote:
>
> Joel,
>
> You should not rely on any particular sorting of the dofs in the solution 
> vector. Instead you can ask the FiniteElement on each cell for the 
> component of a local dof by
> const unsigned int component = fe.system_to_component_index(i).first; 
>
> Apart from that the DataOut object knows about components and you can do 
> some postprocessing on the solution vector using DataPostprocessor [1]
>
> Best,
> Daniel
>
>  [1] https://dealii.org/8.4.1/doxygen/deal.II/classDataPostprocessor.html
>
> Am Montag, 22. August 2016 17:05:08 UTC+2 schrieb Joel Davidsson:
>>
>> Dear Daniel,
>>
>> Thank you for a very good answer, adding the fe_values_scalar and 
>> cell_scalar fixed the problem.
>>
>> I have a follow-up question about the solution I get out, how is the data 
>> organized in the solution vector? Say for example I want to loop over all 
>> the x components, how would I do that?
>>
>> Best,
>> Joel
>>
>> On Friday, August 12, 2016 at 5:35:07 PM UTC+2, Daniel Arndt wrote:
>>>
>>> Joel,
>>>
>>> Yes the matrix you are assembling is a vector-valued mass matrix now.
>>> For me your code is failing with
>>>
>>> void dealii::FEValuesBase::get_function_gradients(const 
>>> InputVector&, std::vector>> InputVector::value_type> >&) const [with InputVector = 
>>> dealii::Vector; int dim = 3; int spacedim = 3; typename 
>>> InputVector::value_type = double]
>>> The violated condition was: 
>>> (fe->n_components()) == (1)
>>> The name and call sequence of the exception was:
>>> dealii::ExcDimensionMismatch((fe->n_components()),(1))
>>> Additional Information: 
>>> Dimension 3 not equal to 1
>>>
>>> and this is not surprising. What you are doing is to extract the 
>>> gradients of a vector-valued finite element solution. The object that 
>>> should store this should therefore be a 
>>> std::vector> as you want to store for each 
>>> quadrature point and each component a Tensor<1,dim>.
>>>
>>> What you want to do is really that in "Do not work". As you have two 
>>> DoFHandlers, you should also have two FEValues objects and two 
>>> cell_iterator corresponding to the correct DoFHandler. Then you can extract 
>>> the gradient of your scalar field and project it onto the ansatz space for 
>>> the vector-valued field. This should look like:
>>>
>>> typename DoFHandler::active_cell_iterator cell_scalar = 
>>> dof_handler_scalar.begin_active();
>>> typename DoFHandler::active_cell_iterator cell = 
>>> dof_handler.begin_active();
>>> ...
>>>
>>>  for (; cell!=endc; ++cell, ++cell_vector)
>>> {
>>>   fe_values.reinit (cell);
>>>   fe_values_scalar.reinit (cell_scalar);
>>>   fe_values_scalar.get_function_gradients(test,fe_function_grad);
>>>   ...
>>> }
>>>
>>> Best,
>>> Daniel
>>>
>>

-- 
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 dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Vector-valued gradient of solution vector

[Daniel Arndt]

Joel,

You should not rely on any particular sorting of the dofs in the solution 
vector. Instead you can ask the FiniteElement on each cell for the 
component of a local dof by
const unsigned int component = fe.system_to_component_index(i).first; 

Apart from that the DataOut object knows about components and you can do 
some postprocessing on the solution vector using DataPostprocessor [1]

Best,
Daniel

 [1] https://dealii.org/8.4.1/doxygen/deal.II/classDataPostprocessor.html

Am Montag, 22. August 2016 17:05:08 UTC+2 schrieb Joel Davidsson:
>
> Dear Daniel,
>
> Thank you for a very good answer, adding the fe_values_scalar and 
> cell_scalar fixed the problem.
>
> I have a follow-up question about the solution I get out, how is the data 
> organized in the solution vector? Say for example I want to loop over all 
> the x components, how would I do that?
>
> Best,
> Joel
>
> On Friday, August 12, 2016 at 5:35:07 PM UTC+2, Daniel Arndt wrote:
>>
>> Joel,
>>
>> Yes the matrix you are assembling is a vector-valued mass matrix now.
>> For me your code is failing with
>>
>> void dealii::FEValuesBase::get_function_gradients(const 
>> InputVector&, std::vector> InputVector::value_type> >&) const [with InputVector = 
>> dealii::Vector; int dim = 3; int spacedim = 3; typename 
>> InputVector::value_type = double]
>> The violated condition was: 
>> (fe->n_components()) == (1)
>> The name and call sequence of the exception was:
>> dealii::ExcDimensionMismatch((fe->n_components()),(1))
>> Additional Information: 
>> Dimension 3 not equal to 1
>>
>> and this is not surprising. What you are doing is to extract the 
>> gradients of a vector-valued finite element solution. The object that 
>> should store this should therefore be a 
>> std::vector> as you want to store for each 
>> quadrature point and each component a Tensor<1,dim>.
>>
>> What you want to do is really that in "Do not work". As you have two 
>> DoFHandlers, you should also have two FEValues objects and two 
>> cell_iterator corresponding to the correct DoFHandler. Then you can extract 
>> the gradient of your scalar field and project it onto the ansatz space for 
>> the vector-valued field. This should look like:
>>
>> typename DoFHandler::active_cell_iterator cell_scalar = 
>> dof_handler_scalar.begin_active();
>> typename DoFHandler::active_cell_iterator cell = 
>> dof_handler.begin_active();
>> ...
>>
>>  for (; cell!=endc; ++cell, ++cell_vector)
>> {
>>   fe_values.reinit (cell);
>>   fe_values_scalar.reinit (cell_scalar);
>>   fe_values_scalar.get_function_gradients(test,fe_function_grad);
>>   ...
>> }
>>
>> Best,
>> Daniel
>>
>

-- 
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 dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Vector-valued gradient of solution vector

[Joel Davidsson]

Dear Daniel,

Thank you for a very good answer, adding the fe_values_scalar and 
cell_scalar fixed the problem.

I have a follow-up question about the solution I get out, how is the data 
organized in the solution vector? Say for example I want to loop over all 
the x components, how would I do that?

Best,
Joel

On Friday, August 12, 2016 at 5:35:07 PM UTC+2, Daniel Arndt wrote:
>
> Joel,
>
> Yes the matrix you are assembling is a vector-valued mass matrix now.
> For me your code is failing with
>
> void dealii::FEValuesBase::get_function_gradients(const 
> InputVector&, std::vector InputVector::value_type> >&) const [with InputVector = 
> dealii::Vector; int dim = 3; int spacedim = 3; typename 
> InputVector::value_type = double]
> The violated condition was: 
> (fe->n_components()) == (1)
> The name and call sequence of the exception was:
> dealii::ExcDimensionMismatch((fe->n_components()),(1))
> Additional Information: 
> Dimension 3 not equal to 1
>
> and this is not surprising. What you are doing is to extract the gradients 
> of a vector-valued finite element solution. The object that should store 
> this should therefore be a std::vector> as you 
> want to store for each quadrature point and each component a Tensor<1,dim>.
>
> What you want to do is really that in "Do not work". As you have two 
> DoFHandlers, you should also have two FEValues objects and two 
> cell_iterator corresponding to the correct DoFHandler. Then you can extract 
> the gradient of your scalar field and project it onto the ansatz space for 
> the vector-valued field. This should look like:
>
> typename DoFHandler::active_cell_iterator cell_scalar = 
> dof_handler_scalar.begin_active();
> typename DoFHandler::active_cell_iterator cell = 
> dof_handler.begin_active();
> ...
>
>  for (; cell!=endc; ++cell, ++cell_vector)
> {
>   fe_values.reinit (cell);
>   fe_values_scalar.reinit (cell_scalar);
>   fe_values_scalar.get_function_gradients(test,fe_function_grad);
>   ...
> }
>
> Best,
> Daniel
>

-- 
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 dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Vector-valued gradient of solution vector

[Daniel Arndt]

Joel,

Yes the matrix you are assembling is a vector-valued mass matrix now.
For me your code is failing with

void dealii::FEValuesBase::get_function_gradients(const 
InputVector&, std::vector >&) const [with InputVector = 
dealii::Vector; int dim = 3; int spacedim = 3; typename 
InputVector::value_type = double]
The violated condition was: 
(fe->n_components()) == (1)
The name and call sequence of the exception was:
dealii::ExcDimensionMismatch((fe->n_components()),(1))
Additional Information: 
Dimension 3 not equal to 1

and this is not surprising. What you are doing is to extract the gradients 
of a vector-valued finite element solution. The object that should store 
this should therefore be a std::vector> as you 
want to store for each quadrature point and each component a Tensor<1,dim>.

What you want to do is really that in "Do not work". As you have two 
DoFHandlers, you should also have two FEValues objects and two 
cell_iterator corresponding to the correct DoFHandler. Then you can extract 
the gradient of your scalar field and project it onto the ansatz space for 
the vector-valued field. This should look like:

typename DoFHandler::active_cell_iterator cell_scalar = 
dof_handler_scalar.begin_active();
typename DoFHandler::active_cell_iterator cell = 
dof_handler.begin_active();
...

 for (; cell!=endc; ++cell, ++cell_vector)
{
  fe_values.reinit (cell);
  fe_values_scalar.reinit (cell_scalar);
  fe_values_scalar.get_function_gradients(test,fe_function_grad);
  ...
}

Best,
Daniel

-- 
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 dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Vector-valued gradient of solution vector

[Joel Davidsson]

Dear Daniel,

Thank you for a fast answer. I think I got the matrix and right-hand side 
correct now, please see test.cc for the example code.

What I want to from the start was to get the vector-valued gradient from a 
scalar field, but I could not get this to work. See section "Do not work" 
in function setup_system in test.cc, this gives 0 everywhere. I want to use 
the information from test vector but since it was made with a different 
dof_handler it doesn´t seem to work. Is there a way around this problem? 
Can a vector be transferred from one dof_handler to another? or is there a 
better way to solve this?

Best,
Joel

On Friday, August 12, 2016 at 12:00:36 AM UTC+2, Daniel Arndt wrote:
>
> Joel,
>
> Does the equation you want to solve for have multiple components? 
> Otherwise it would not make sense to multiply in the right hand side 
> something with a gradient.
> If you want to project the gradient into a vector valued finite element 
> ansatz space, then the matrix you are assembling looks weird.
> In what you are doing you are also considering the coupling between all 
> components.
>
> You are calculating the gradient of your scalar finite element function 
> correctly, but you need to find the correct component to use.
> For doing this you can use something like:
>
> const unsigned int component = fe.system_to_component_index(i).first;
> cell_rhs(i) += ((fe_values.shape_value (i, q_index) *
> fe_function_grad[q_index][component])*
> fe_values.JxW (q_index));
>
> If you want to assemble a mass matrix for a vector-valued finite element, 
> you have also to restrict assembling for the matrix to the case
> where the component for the ith and jth ansatz function are the same.
>
> You might want to have a look at step-8 [1] and the module "Handling 
> vector valued problems"[2].
>
> Best,
> Daniel
>
> [1] 
> https://www.dealii.org/8.4.1/doxygen/deal.II/step_8.html#ElasticProblemassemble_system
> [2] 
> https://www.dealii.org/8.4.1/doxygen/deal.II/group__vector__valued.html
>

-- 
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 dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
/* -
 *
 * Copyright (C) 1999 - 2013 by the deal.II authors
 *
 * This file is part of the deal.II library.
 *
 * The deal.II library is free software; you can use it, redistribute
 * it, and/or modify it under the terms of the GNU Lesser General
 * Public License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 * The full text of the license can be found in the file LICENSE at
 * the top level of the deal.II distribution.
 *
 * -

 *
 * Author: Wolfgang Bangerth, University of Heidelberg, 1999
 */

#include 

#include 
#include 
#include 

#include 
#include 
#include 

#include 
#include 
#include 
#include 

#include 
#include 

#include 
#include 

#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 

#include 
#include 
#include 

#include 

using namespace dealii;

const double pi=3.14159265359;
double center_id=0;

template 
class Step4
{
public:
  Step4 ();
  void run ();

private:
  void make_grid ();
  void setup_system();
  void assemble_system ();
  void solve ();
  void output_results () const;

  Triangulation   triangulation;
  FE_Qfe_scalar;
  FESystem		   fe;

  DoFHandler  dof_handler_scalar;
  DoFHandler  dof_handler;

  SparsityPattern  sparsity_pattern;
  SparseMatrix system_matrix;

  ConstraintMatrix constraints;

  Vector   solution;
  Vector   system_rhs;
  Vector   test;

  std::vector> nodes;
};

template 
Step4::Step4 ()
  :
  fe_scalar (1),
  fe (FE_Q(1),dim),
  dof_handler (triangulation),
  dof_handler_scalar (triangulation)
{}

template 
void Step4::make_grid ()
{
  GridGenerator::hyper_cube (triangulation, -2.5, 2.5);
  triangulation.refine_global (4);

  for (unsigned int step=0; step<2; ++step)
  {
	  typename Triangulation::active_cell_iterator
	  cell = triangulation.begin_active(),
	  endc = triangulation.end();
	  for (; cell!=endc; ++cell)
	  {
		  for (unsigned int v=0;v < GeometryInfo::vertices_per_cell;++v)
		  {
			  Point center;
			  for (unsigned int i=0; ivertex(v));
			  if (std::fabs(distance_from_center) < 1.0)
			  {
  cell->set_refine_flag ();
  break;
			  }
		  }
	  }
	  triangulation.execute_coarsening_and_refinement ();
  }

  std::cout << "   Number of active cells: "
<< triangulation.n_active_cells()
<< std::endl
<< "   Total 

[deal.II] Re: Vector-valued gradient of solution vector

[Daniel Arndt]

Joel,

Does the equation you want to solve for have multiple components? Otherwise 
it would not make sense to multiply in the right hand side something with a 
gradient.
If you want to project the gradient into a vector valued finite element 
ansatz space, then the matrix you are assembling looks weird.
In what you are doing you are also considering the coupling between all 
components.

You are calculating the gradient of your scalar finite element function 
correctly, but you need to find the correct component to use.
For doing this you can use something like:

const unsigned int component = fe.system_to_component_index(i).first;
cell_rhs(i) += ((fe_values.shape_value (i, q_index) *
fe_function_grad[q_index][component])*
fe_values.JxW (q_index));

If you want to assemble a mass matrix for a vector-valued finite element, 
you have also to restrict assembling for the matrix to the case
where the component for the ith and jth ansatz function are the same.

You might want to have a look at step-8 [1] and the module "Handling vector 
valued problems"[2].

Best,
Daniel

[1] 
https://www.dealii.org/8.4.1/doxygen/deal.II/step_8.html#ElasticProblemassemble_system
[2] https://www.dealii.org/8.4.1/doxygen/deal.II/group__vector__valued.html

-- 
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 dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.