Re: [deal.II] Re: How to get the coordinates of a given degree of freedom on an edge?

2019-02-18 Thread Phạm Ngọc Kiên
Dear Prof. Wolfgang Bangerth,
I would like to thank you very much.
I am reading through the FE_NedelecSZ codes in Deal.II library.
I currently try to solve the electromagnetics problem, followed by the
research of Ross Kynch.
My goal is to compute the right hand side for the problem:
curl curl E + i*omega*mu*sigma*E = - Js
I only want to set the Js term here to be none zero, for example to be 1,
at a given degree of freedom (i.e. at source position), and  0 everywhere
else.
Thus, I need to write the codes to query if a degree of freedom is at
source position.
As you mentioned, now I think that I need to check if the source point is
on an edge when assembling the right hand side vector.
Could you please tell me help me about this?
Best regards,
Pham Ngoc Kien

Vào Th 3, 19 thg 2, 2019 vào lúc 15:45 Wolfgang Bangerth <
bange...@colostate.edu> đã viết:

> On 2/18/19 10:47 PM, Phạm Ngọc Kiên wrote:
> >
> > 2. I have tried some thing with the third method. And here-below is my
> code:
> >
> > Point p{0.5,0.5,0.5};//position in reference cell
> > Quadrature q(p);
> > FEValues fe_values_q(fe, q,update_quadrature_points);
> > fe_values_q.reinit(cell);
> >
> > //position in the real cell
> > std::vector> dof_pos = fe_values_q.get_quadrature_points();
> > for (int k =0; k  >  for (int j =0; j  >  std::cout< >  std::cout<<"\n";
> >  }
> > }
> >
> >  From this codes, I can get the real position of arbitrary points in my
> current reference cell.
> >
> > However, my purpose is to get the real position of a given degree of
> freedom.
> >
> > Could you please tell me how to get the position of a degree of freedom
> in the reference cell?
>
> The problem is that the NedelecSZ element does not have any such position.
> You
> are thinking of Lagrange elements (or the regular Nedelec elements) for
> which
> every shape function has a specific location where it is defined: It is
> one
> there, and zero at all other such locations. This is what we call "support
> points" in deal.II.
>
> But the NedelecSZ element is not defined this way (at least that's what I
> believe to be the case). Rather, shape functions are defined in such a way
> that certain *integrals* (rather than *point values*) are zero or one for
> all
> shape functions. As a consequence, there are no support points: Shape
> functions (and degrees of freedom) are not defined at individual points,
> but
> at whole edges.
>
> This means that your approach can not work. I don't know what you want to
> do,
> but let's start with that: What is your *goal*? Maybe then we can think
> about
> how to achieve it!
>
> Best
>   W.
>
> --
> 
> Wolfgang Bangerth  email: bange...@colostate.edu
> 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 dealii+unsubscr...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
>

-- 
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.


Re: [deal.II] Re: How to get the coordinates of a given degree of freedom on an edge?

2019-02-18 Thread Wolfgang Bangerth
On 2/18/19 10:47 PM, Phạm Ngọc Kiên wrote:
> 
> 2. I have tried some thing with the third method. And here-below is my code:
> 
> Point p{0.5,0.5,0.5};//position in reference cell
> Quadrature q(p);
> FEValues fe_values_q(fe, q,update_quadrature_points);
> fe_values_q.reinit(cell);
> 
> //position in the real cell
> std::vector> dof_pos = fe_values_q.get_quadrature_points();
> for (int k =0; k   for (int j =0; j   std::cout<  std::cout<<"\n";
>  }
> }
> 
>  From this codes, I can get the real position of arbitrary points in my 
> current reference cell.
> 
> However, my purpose is to get the real position of a given degree of freedom.
> 
> Could you please tell me how to get the position of a degree of freedom in 
> the reference cell?

The problem is that the NedelecSZ element does not have any such position. You 
are thinking of Lagrange elements (or the regular Nedelec elements) for which 
every shape function has a specific location where it is defined: It is one 
there, and zero at all other such locations. This is what we call "support 
points" in deal.II.

But the NedelecSZ element is not defined this way (at least that's what I 
believe to be the case). Rather, shape functions are defined in such a way 
that certain *integrals* (rather than *point values*) are zero or one for all 
shape functions. As a consequence, there are no support points: Shape 
functions (and degrees of freedom) are not defined at individual points, but 
at whole edges.

This means that your approach can not work. I don't know what you want to do, 
but let's start with that: What is your *goal*? Maybe then we can think about 
how to achieve it!

Best
  W.

-- 

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


Re: [deal.II] Re: How to get the coordinates of a given degree of freedom on an edge?

2019-02-18 Thread Phạm Ngọc Kiên
I am sorry as my question is not clear enough.
What I mean is when we loop over all degree of freedom in a cell, for
example,

for (unsigned int i = 0; i < dofs_per_cell; ++i) { .. }

Do we have any method to get the position of the i_th degree of
freedom inside this loop?

I would like to thank you very much!

Best,

Pham Ngoc Kien


Vào Th 3, 19 thg 2, 2019 vào lúc 14:47 Phạm Ngọc Kiên <
ngockien.lea...@gmail.com> đã viết:

> Hi,
> I have tried with the methods described in the FAQ and figured out that:
> 1. As I used edge-based element (FE_NedelecSZ), I queried if my element
> had support points by  *fe.has_support_points()*. As a consequence, it
> returned zero.
> The *FiniteElement::get_unit_support_points() *also returned an empty
> vector.
> As there was no support point, DoFTools::map_dofs_to_support_points()
> could not work as you mentioned above.
>
> 2. I have tried some thing with the third method. And here-below is my
> code:
>
> Point p{0.5,0.5,0.5}; //position in reference cell
> Quadrature q(p);
> FEValues fe_values_q(fe, q, update_quadrature_points);
> fe_values_q.reinit(cell);
>
> //position in the real cell
> std::vector> dof_pos = fe_values_q.get_quadrature_points();
> for (int k = 0; k  for (int j = 0; j  std::cout< std::cout<<"\n";
> }
> }
>
> From this codes, I can get the real position of arbitrary points in my 
> current reference cell.
>
> However, my purpose is to get the real position of a given degree of freedom.
>
> Could you please tell me how to get the position of a degree of freedom in 
> the reference cell?
>
> Thank you very much for your help.
>
> Best,
>
> Pham Ngoc Kien
>
>
> Vào Th 3, 19 thg 2, 2019 vào lúc 11:47 Phạm Ngọc Kiên <
> ngockien.lea...@gmail.com> đã viết:
>
>> Yes, thank you very much.
>> I will dig deeper into this.
>> Best,
>> Pham Ngoc Kien
>>
>> Vào Th 3, 19 thg 2, 2019 vào lúc 11:23 Bruno Turcksin <
>> bruno.turck...@gmail.com> đã viết:
>>
>>> Pham Ngoc Kien,
>>>
>>> Le lun. 18 févr. 2019 à 19:33, Phạm Ngọc Kiên 
>>> a écrit :
>>>
 With the DoFTools::map_dofs_to_support_points() , I saw in the library
 manual that it can be used with no edge elements or the like.
 As my finite element is edge-based one, I wonder if this function work?

>>> No the function won't work but there are two other methods described in
>>> the FAQ that you could use. In particular, you can create a special
>>> quadrature at the points of interest.
>>>
>>> Best,
>>>
>>> Bruno
>>>
>>> --
>>> 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.
>>>
>>

-- 
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.


Re: [deal.II] Re: How to get the coordinates of a given degree of freedom on an edge?

2019-02-18 Thread Phạm Ngọc Kiên
Hi,
I have tried with the methods described in the FAQ and figured out that:
1. As I used edge-based element (FE_NedelecSZ), I queried if my element had
support points by  *fe.has_support_points()*. As a consequence, it returned
zero.
The *FiniteElement::get_unit_support_points() *also returned an empty
vector.
As there was no support point, DoFTools::map_dofs_to_support_points() could
not work as you mentioned above.

2. I have tried some thing with the third method. And here-below is my code:

Point p{0.5,0.5,0.5}; //position in reference cell
Quadrature q(p);
FEValues fe_values_q(fe, q, update_quadrature_points);
fe_values_q.reinit(cell);

//position in the real cell
std::vector> dof_pos = fe_values_q.get_quadrature_points();
for (int k = 0; k From this codes, I can get the real position of arbitrary points in my
current reference cell.

However, my purpose is to get the real position of a given degree of freedom.

Could you please tell me how to get the position of a degree of
freedom in the reference cell?

Thank you very much for your help.

Best,

Pham Ngoc Kien


Vào Th 3, 19 thg 2, 2019 vào lúc 11:47 Phạm Ngọc Kiên <
ngockien.lea...@gmail.com> đã viết:

> Yes, thank you very much.
> I will dig deeper into this.
> Best,
> Pham Ngoc Kien
>
> Vào Th 3, 19 thg 2, 2019 vào lúc 11:23 Bruno Turcksin <
> bruno.turck...@gmail.com> đã viết:
>
>> Pham Ngoc Kien,
>>
>> Le lun. 18 févr. 2019 à 19:33, Phạm Ngọc Kiên 
>> a écrit :
>>
>>> With the DoFTools::map_dofs_to_support_points() , I saw in the library
>>> manual that it can be used with no edge elements or the like.
>>> As my finite element is edge-based one, I wonder if this function work?
>>>
>> No the function won't work but there are two other methods described in
>> the FAQ that you could use. In particular, you can create a special
>> quadrature at the points of interest.
>>
>> Best,
>>
>> Bruno
>>
>> --
>> 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.
>>
>

-- 
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.


Re: [deal.II] Re: How to get the coordinates of a given degree of freedom on an edge?

2019-02-18 Thread Phạm Ngọc Kiên
Yes, thank you very much.
I will dig deeper into this.
Best,
Pham Ngoc Kien

Vào Th 3, 19 thg 2, 2019 vào lúc 11:23 Bruno Turcksin <
bruno.turck...@gmail.com> đã viết:

> Pham Ngoc Kien,
>
> Le lun. 18 févr. 2019 à 19:33, Phạm Ngọc Kiên 
> a écrit :
>
>> With the DoFTools::map_dofs_to_support_points() , I saw in the library
>> manual that it can be used with no edge elements or the like.
>> As my finite element is edge-based one, I wonder if this function work?
>>
> No the function won't work but there are two other methods described in
> the FAQ that you could use. In particular, you can create a special
> quadrature at the points of interest.
>
> Best,
>
> Bruno
>
> --
> 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.
>

-- 
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.


Re: [deal.II] Re: How to get the coordinates of a given degree of freedom on an edge?

2019-02-18 Thread Bruno Turcksin
Pham Ngoc Kien,

Le lun. 18 févr. 2019 à 19:33, Phạm Ngọc Kiên  a
écrit :

> With the DoFTools::map_dofs_to_support_points() , I saw in the library
> manual that it can be used with no edge elements or the like.
> As my finite element is edge-based one, I wonder if this function work?
>
No the function won't work but there are two other methods described in the
FAQ that you could use. In particular, you can create a special quadrature
at the points of interest.

Best,

Bruno

-- 
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.


Re: [deal.II] Re: How to get the coordinates of a given degree of freedom on an edge?

2019-02-18 Thread Phạm Ngọc Kiên
Hi Bruno,
With the DoFTools::map_dofs_to_support_points() , I saw in the library
manual that it can be used with no edge elements or the like.
As my finite element is edge-based one, I wonder if this function work?
I have check the function fe.has_support_points() in my code and it
returned 0.
That is the reason why I put my question.
Thank you very much for your suggestion.

Best,
Pham Ngoc Kien

Vào Th 3, 19 thg 2, 2019 vào lúc 06:18 Bruno Turcksin <
bruno.turck...@gmail.com> đã viết:

> Hi,
>
> Take a look at this entry
> 
> in the FAQ.
>
> Best,
>
> Bruno
>
> On Monday, February 18, 2019 at 12:20:48 AM UTC-5, Phạm Ngọc Kiên wrote:
>>
>> Dear colleagues,
>> I am working with FE_Nedelec_SZ element (edge-based finite element) in
>> Deal.II (version 9.0.0).
>> For my current problem, I need to get the global coordinates of a given
>> degree of freedom on and edge.
>> Do we have any code for this task?
>>
>> I have also thought of using material_id for an edge as an alternative
>> way to address my problem. However, I personally think that the material_id
>> can be used only with a cell.
>>
>> I would like to thank you very much in advance.
>>
>> Pham Ngoc Kien
>> Ph.D. student
>> Seoul National University
>>
> --
> 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.
>

-- 
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: Sparsematrix initialization in P4est program

2019-02-18 Thread David Wells
Hi Gabriel,

I believe that you are still using the standard deal.II SparseMatrix class; 
is this correct? If so, this class won't work with fully distributed 
calculations. You will need to use either the Trilinos or PETSc wrappers 
since not all information is available on the local processor with the 
default SparseMatrix class.

Either way: I believe the error message you encountered is a result of a 
bug in deal.II. I don't think that the routine 
GridTools::get_finest_common_cells(), which is called by that particular 
make_sparsity_pattern variant, works with distributed triangulations. If 
you want to use distributed triangulations you will need to set up the 
coupled sparsity pattern in a different way.

If you want to parallelize step-35, it may be better to start by looking at 
step-32, which shows how to implement the Stokes equations (coupled to a 
Boussinesq equation) in a fully distributed setting. I believe they set up 
the coupling of the matrix with velocity and pressure blocks in a different 
way that should work for you.

Does this make sense?

Thanks,
David Wells



On Sunday, February 17, 2019 at 4:29:06 PM UTC-5, gabriel...@koeln.de wrote:
>
> Hey everybody, 
> I have a problem with initializing a Sparsematrix in a parallel program.
>  am using a parallel::distributed::triangulaton (as in the Step40 
> tutorial) on 
> the Navier Stokes Projection solver from the step35 tutorial.
> The problem arises in the function 
>
> NavierStokesProjection::initialize_gradient_operator()
> {
> {
> DynamicSparsityPattern 
> 
>  
> dsp(dof_handler_velocity.n_dofs(), dof_handler_pressure.n_dofs());
> DoFTools::make_sparsity_pattern 
> 
>  
> (dof_handler_velocity, dof_handler_pressure, dsp);
> sparsity_pattern_pres_vel.copy_from (dsp);
> }
>
> where I have to use two different FE-spaces, one for the pressure and one 
> for the velocity.
> To perform this in parallel I used the following code (according to the 
> sparsematrices in step40)
>
> {
> DynamicSparsityPattern 
> dsp(locally_relevant_dofs_vel.size(),locally_relevant_dofs_pres.size());
> DoFTools::make_sparsity_pattern (dof_handler_velocity, 
> dof_handler_pressure, dsp); 
> sparsity_pattern_pres_vel.copy_from (dsp);
> }
>
> This code compilates fine, but when Running it with MPI, it stops in the 
> "make_sparsity_pattern" line with the error code
> "this->is_artificial() == false
> Additional information:
> Cant asks for Dofs on articial cells."
>
> Does somebody know how I can fix this? 
>
> My last try was to use the Indexsets instead of their size, when declaring 
> dsp.
> But unfortunately DynamicSparsityPattern can only take on IndexSet.
>
>
> Thanks a lot
>
> Gabriel
>
>
>
>
>

-- 
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: How to get the coordinates of a given degree of freedom on an edge?

2019-02-18 Thread Bruno Turcksin
Hi,

Take a look at this entry 

 
in the FAQ.

Best,

Bruno

On Monday, February 18, 2019 at 12:20:48 AM UTC-5, Phạm Ngọc Kiên wrote:
>
> Dear colleagues,
> I am working with FE_Nedelec_SZ element (edge-based finite element) in 
> Deal.II (version 9.0.0).
> For my current problem, I need to get the global coordinates of a given 
> degree of freedom on and edge. 
> Do we have any code for this task?
>
> I have also thought of using material_id for an edge as an alternative way 
> to address my problem. However, I personally think that the material_id can 
> be used only with a cell.
>
> I would like to thank you very much in advance.
>
> Pham Ngoc Kien
> Ph.D. student
> Seoul National University
>

-- 
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: Calculate error values as percentages

2019-02-18 Thread Daniel Arndt
Maxi,

Usually when calculating error values I use 
> VectorTools::integrate_difference(). Still, that gives me an absolute 
> value. Now I would like to calculate the relative value, i.e. 
> |u-u_h|_{H^1}/|u|_{H^1}. Based on what I could find until now that means 
> that I have to create a vector containing the value for the gradients of 
> the correct solution at each point (using 
> VectorTools::get_position_vector() and a loop afterwards), which then can 
> be used to calculate |u|_{H^1} by calling vec.h1_norm(). The problem is, 
> though, that evaluating the solution at all points could be quite costly. 
> Thus I was wondering if there are any other functions which could help 
> speed it up?
>
 
In general, you can not calculate |u|_{H^1} exactly at all. All you can do, 
is to approximate the integrl using some quadrature rule.
In particular, the norm of the continuous solution is independent of the 
discretization used. Hence, you can use an entirely separate mesh (or just 
provide the norm from the outside).
Normally, calculating the relative error with respect to the discrete 
solution (|u-u_h|_{H^1}/|u_h|_{H^1}) should also suffice since you expect 
the error to be small anyway.

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] Calculate error values as percentages

2019-02-18 Thread 'Maxi Miller' via deal.II User Group
Usually when calculating error values I use 
VectorTools::integrate_difference(). Still, that gives me an absolute 
value. Now I would like to calculate the relative value, i.e. 
|u-u_h|_{H^1}/|u|_{H^1}. Based on what I could find until now that means 
that I have to create a vector containing the value for the gradients of 
the correct solution at each point (using 
VectorTools::get_position_vector() and a loop afterwards), which then can 
be used to calculate |u|_{H^1} by calling vec.h1_norm(). The problem is, 
though, that evaluating the solution at all points could be quite costly. 
Thus I was wondering if there are any other functions which could help 
speed it up?

-- 
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.