Re: [deal.II] Finding the y-coordinate of the boundary during the assembly

[Wolfgang Bangerth]

On 3/14/19 5:50 PM, jane@jandj-ltd.com wrote:
> 
> I've tried looking a cell->face(face_no)->at_boundary() type input within the 
> cell iterator and by using something like cell->vertex(v)(1) but I can't 
> quite 
> get my head around how to find the y coordinate of the top boundary for every 
> cell/quadrature point calculation.

Jane,
if you've already identified that a face is at the boundary, then
   cell->face(face_no)->center()
returns a Point object. So if you need the y-coordinate, it would be
   cell->face(face_no)->center()[1]

Or are you asking the following question: "For a given quadrature point at 
(x,y), find how far the domain extends above (x,y) in y-direction"?

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] Accuracy of Dirichlet condition for p in step-20

[Wolfgang Bangerth]

On 3/14/19 5:02 PM, jane@jandj-ltd.com wrote:
> 
> I think just a comment in the bit on how to implement the dirichlet bc in the 
> weak form would be sufficient - something to say 'In the case of an 
> inhomogeneous boundary condition, you would need to set local_rhs = 0 before 
> adding the cell contributions for the boundary condition'.

What line of the code would that be? Do you think it would be wrong to just 
*always* set local_rhs=0?


> I'm still unsure about the step-22 condition, for the normal component for 
> the 
> normal stress. Is this equivalent to a dirichlet condition on the pressure 
> only? I'm a little confused on this one and any thoughts would be helpful.

The normal stress shows up in the weak form after integrating by parts both 
the div(2*eta*eps(u)) and the grad p terms. The term is going to be something 
like
   (v, (2*eta*eps(u) - pI).n)_{Gamma}

So if you prescribe
* no normal flux
* no tangential stress
(i.e., "free slip"), then you will have:
* the normal component of v is zero
* the tangential component of (2*eta*eps(u) - pI).n is zero
As a consequence, the entire boundary term disappears, as you can see if you 
write the product of test function and normal stress as
   (v, s.n) = (v_n, (s.n).n) + (v_t, (s.n) x n)
(normal component plus tangential component of vectors), where the first of 
the two terms on the right disappears because v_n=0 and the second because 
(s.n) x n = 0.

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] ILU preconditioner and mesh-scale dependency?

[Jaekwang Kim]

Thanks for sharing insight! 


On Thursday, March 14, 2019 at 12:18:30 PM UTC-5, Wolfgang Bangerth wrote:
>
> On 3/13/19 8:18 PM, Jaekwang Kim wrote: 
> > 
> > 
> > I wonder now ... 
> > 
> > 1. why smaller scale mesh only results such error, 
>
> You have two terms in your matrix, one from the advection term and one 
> from the reaction term. If you scale the mesh by a factor of c, the 
> height of the shape functions used in the reaction term does not change, 
> but the gradient of shape functions used in the advection term is scaled 
> by 1/c. As a consequence, the balance of the two terms that together 
> make up the matrix changes and you apparently get into a situation where 
> the resulting matrix is ill-conditioned or indeed not invertible. 
>
>
> > 2. If possible, is there any choice of better preconditionner for my 
> > system other than SparseILU to go around this problem? 
>
> To make sure that the system you have assembled is correct, I usually 
> first check with the SparseDirectUMFPACK solver. Once I know that the 
> system is correct, I think about iterative solvers and preconditioners. 
> In your case, the system you have is advection dominated, for which it 
> is notoriously difficult to construct good preconditioners. 
>
> Best 
>   WB 
>
> -- 
>  
> Wolfgang Bangerth  email: bang...@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] what is the different between VectorTools::point_value and fe_values.get_function_values()?

[Phạm Ngọc Kiên]

I highly appreciate all of your answers. It becomes clearly for me to look
through the library.
Let me show all of my function and my ideas below.
I set the point for evaluating its solution to be on the vertex of a cell,
and using the following function to get the solution

template
void OutputTools::point_solution(const Mapping &mapping,
//input mapping
  const DoFHandler
&dof_handler,   //input dof_handler
  const Vector &solution,
 //input solution vector
  const Point &point,
 //input point in real coordinates to get the solution
  Vector &point_solution)
const {  //output

const std::pair::active_cell_iterator, Point>
cell_point =
GridTools::find_active_cell_around_point(mapping, dof_handler, point);
 //find active cell around input point

const double delta = 1e-20;
Assert(GeometryInfo::distance_to_unit_cell(cell_point.second) < delta,
   ExcInternalError());

const Quadrature
quadrature(GeometryInfo::project_to_unit_cell(cell_point.second));
//get quadrature by project the point to unit cell
FEValues fe_values(mapping, dof_handler.get_fe(), quadrature,
update_values | update_gradients |
update_quadrature_points);
fe_values.reinit(cell_point.first);

const FEValuesExtractors::Vector vector_re(0);
const FEValuesExtractors::Vector vector_im(dim);

std::vector> temp_solution_re(quadrature.size());
std::vector> temp_solution_im(quadrature.size());

//Get e_field at point
fe_values[vector_re].get_function_values(solution, temp_solution_re);
fe_values[vector_im].get_function_values(solution, temp_solution_im);
for (unsigned int j = 0; j < dim; ++j) {
point_solution(j) = temp_solution_re[0][j];
point_solution(j + dim) = temp_solution_im[0][j];
}

   //using point_value() for comparing the result

   Vector values (6);
   VectorTools::point_value(dof_handler,solution,point,values);

//print out all the results for testing the code

   std::cout<<"\n";
   std::cout<<"cell: \t"<::active_cell_iterator, Point>
cell_point =
GridTools::find_active_cell_around_point(mapping, dof_handler, point);

I personally think that as I set it to be const, so the return cell
(cell_point.first) is only one cell, rather than a group of cells.
Could you please tell me about this?


2. My second idea is with the codes:
const Quadrature
quadrature(GeometryInfo::project_to_unit_cell(cell_point.second));
FEValues fe_values(mapping, dof_handler.get_fe(), quadrature,
update_values | update_gradients |
update_quadrature_points);
fe_values.reinit(cell_point.first);

I have only one quadrature defined on the cell, and the printing
result show me the same thing.
Can I use this way to get the solution, or it would be more efficient
when using codes like:
QGaussLobatto quadrature_formula(2);
FEValues fe_values(fe, quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
fe_values.reinit(cell_point.first);

And then finding the quadrature point coincide with the evaluating
point to get the solution?

3. My third thought is in my function I have only 1 quadrature point
so that temp_solution_re[0] and temp_solution_im[0] are the variables
containing the solution.

The printing results showed that the solution from point_value() is
identical with those from temp_solution_re[0] and temp_solution_im[0].

Do temp_solution_re[1], temp_solution_re[2],
temp_solution_im[1],temp_solution_im[2] make any sense for my
solution? Because I saw the code can run with these variables.

Thank you very much!


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

> On 3/14/19 1:28 AM, Phạm Ngọc Kiên wrote:
> >
> > So for the first function
> > VectorTools::point_value(dof_handler,solution,point,values);
> > I can specify the observed point in real coordiante, solution vector and
> > dof_handler to get the values, right?
>
> Correct.
>
>
> > For the second one, assuming that I am using QGauss for solving the
> > system, does it mean that if the point is on the edge of the cell then I
> > cannot get the solution at that point?
>
> Yes. It can be any quadrature formula defined on the reference cell. It
> doesn't matter where the points are located.
>
> 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

[deal.II] Finding the y-coordinate of the boundary during the assembly

[jane . lee]

Hi all,

I am solving a problem where one of the functions in the assembly process 
requires the knowledge of the boundary coordinate of the mesh, which moves 
at each timestep. 

My solution should be only varying in the y-direction (it is indeed with 
the current output) and I have a known linear function independent of x so 
I only need the y coordinate of the top boundary above the relevant x-y 
coordinate at each assembly. (I basically need the hydrostatic pressure 
function that is dependent on the depth only).

ie. I need to calculate p = 1.5*(top-boundary-coordinate directly above the 
actually coordinate - coordinate in question) + 1 for each cell vertex in 
the assembly. 


I've tried looking a cell->face(face_no)->at_boundary() type input within 
the cell iterator and by using something like cell->vertex(v)(1) but I 
can't quite get my head around how to find the y coordinate of the top 
boundary for every cell/quadrature point calculation. 

I'd appreciate any help. 

Thanks

-- 
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] Accuracy of Dirichlet condition for p in step-20

[jane . lee]

For step-20,

I think just a comment in the bit on how to implement the dirichlet bc in 
the weak form would be sufficient - something to say 'In the case of an 
inhomogeneous boundary condition, you would need to set local_rhs = 0 
before adding the cell contributions for the boundary condition'.

I'm still unsure about the step-22 condition, for the normal component for 
the normal stress. Is this equivalent to a dirichlet condition on the 
pressure only? I'm a little confused on this one and any thoughts would be 
helpful. thanks!

On Thursday, March 14, 2019 at 10:24:31 PM UTC, Wolfgang Bangerth wrote:
>
>
> Hi Jane, 
>
> > I would think that it might be useful in the tutorial!! 
>
> The same applies here: I deal with too much email to be clear about what 
> "it" replies to :-) If you think "it" would be useful, can you make a 
> concrete suggestion? 
>
>
> > Please do let me know if you have any suggestions on a perhaps analogous 
> > question for step-22 in my other reply to this thread. I'm making the 
> > analogy between the dirichlet boundary condition implementation and the 
> > normal component of the normal stress condition for the partial boundary 
> > condition in step-22. 
>
> The same would apply there as well, if it isn't already done there. So 
> if you think that something ought to be changed in step-22, then I'll be 
> happy to get *concrete* suggestions for that as well! 
>
> Thanks 
>   W. 
>
> -- 
>  
> Wolfgang Bangerth  email: bang...@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] Accuracy of Dirichlet condition for p in step-20

[Wolfgang Bangerth]

Hi Jane,

> I would think that it might be useful in the tutorial!!

The same applies here: I deal with too much email to be clear about what 
"it" replies to :-) If you think "it" would be useful, can you make a 
concrete suggestion?


> Please do let me know if you have any suggestions on a perhaps analogous 
> question for step-22 in my other reply to this thread. I'm making the 
> analogy between the dirichlet boundary condition implementation and the 
> normal component of the normal stress condition for the partial boundary 
> condition in step-22.

The same would apply there as well, if it isn't already done there. So 
if you think that something ought to be changed in step-22, then I'll be 
happy to get *concrete* suggestions for that as well!

Thanks
  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] Accuracy of Dirichlet condition for p in step-20

[Wolfgang Bangerth]

On 3/14/19 11:42 AM, jane@jandj-ltd.com wrote:
> 
> as an addition to this post, For Dirichlet condition that are in the 
> weak form, would this have to be done?
> 
> And for example, for an inhomogeneous normal component of the normal 
> stress condition in step-22 for Stokes (1st of the partial bc discussion 
> in the tutorial), do you also have to initialise local_rhs? I have a 
> similar condition for the top boundary with the Stokes equation but 
> can't quite figure out if the pressure should equal the normal component 
> of the normal stress or not (if the stress tensor/velocity is non-zero, 
> then it shouldn't equal the inhomogeneous condition input?)

I'm not entirely sure I know what you mean with "it" or "this" -- I 
suspect you are asking whether local_rhs should be set to zero every 
time we move to a new cell.

The answer is yes. We think of these as each cell's contribution to the 
global matrix. We want to compute it separately on each cell, rather 
than carry over some state from the previous cell.

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] Accuracy of Dirichlet condition for p in step-20

[jane . lee]

Hi Wolfgang,

I would think that it might be useful in the tutorial!!

Please do let me know if you have any suggestions on a perhaps analogous 
question for step-22 in my other reply to this thread. I'm making the 
analogy between the dirichlet boundary condition implementation and the 
normal component of the normal stress condition for the partial boundary 
condition in step-22. 

Thank you!


On Thursday, March 14, 2019 at 5:12:50 PM UTC, Wolfgang Bangerth wrote:
>
> On 3/14/19 9:13 AM, jane...@jandj-ltd.com  wrote: 
> > Note that I think this was the issue. 
> > 
> > You would need to put local_rhs = 0 within the part below to apply the 
> > Dirichlet condition. Please do let me know if you think this is wrong 
> > but I am getting the correct output. Thanks 
>
> That seems reasonable and correct. 
>
> Are you suggesting that step-20 could be made clearer by adding a line 
> that in that program doesn't make a difference but would if it were used 
> in a more general context? If so, could you suggest what to put and where? 
>
> Best & glad to hear that you solved the problem! 
>   W. 
>
> -- 
>  
> Wolfgang Bangerth  email: bang...@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] Accuracy of Dirichlet condition for p in step-20

[jane . lee]

Hi all,

as an addition to this post, For Dirichlet condition that are in the weak 
form, would this have to be done?

And for example, for an inhomogeneous normal component of the normal stress 
condition in step-22 for Stokes (1st of the partial bc discussion in the 
tutorial), do you also have to initialise local_rhs? I have a similar 
condition for the top boundary with the Stokes equation but can't quite 
figure out if the pressure should equal the normal component of the normal 
stress or not (if the stress tensor/velocity is non-zero, then it shouldn't 
equal the inhomogeneous condition input?)

Thanks


On Thursday, March 14, 2019 at 3:13:59 PM UTC, jane...@jandj-ltd.com wrote:
>
> Note that I think this was the issue. 
>
> You would need to put local_rhs = 0 within the part below to apply the 
> Dirichlet condition. Please do let me know if you think this is wrong but I 
> am getting the correct output. Thanks
>
> for (unsigned int face_no=0;
> face_no::faces_per_cell;
> ++face_no)
> if (cell->face(face_no)->at_boundary()
> &&
> (cell->face(face_no)->boundary_id() == 1)) // top has boundary id 1 
> On Thursday, March 14, 2019 at 2:24:55 PM UTC, jane...@jandj-ltd.com 
> wrote:
>>
>> Hi Wolfgang, Sorry for the reply - for some reason i didn't get the 
>> notification. 
>>
>> What i mean by correct is that I do have two uncoupled equations, but 
>> that it is correct in that I had already verified it using what you 
>> suggested (by using exact solutions - MMS).
>>
>> By it fails, I mean that the top pf value output is not equal to pr as it 
>> should be, which I set in the weak form as in my first message in my 
>> thread: 
>>
>> local_rhs(i) += -( fe_face_values_pf[velocities].value (i, q) *
>> fe_face_values_pf.normal_vector(q) *
>> pr_boundary_values[q] *
>> fe_face_values_pf.JxW(q));
>>
>>  this is in the same way as step-20.
>>
>> I have now done further tests by just letting pr_boundary_values[q] to 
>> 1.0 in the blurb copied here - I do not get 1.0 as my value for the 
>> solution for pf  on boundary with id 1, which I believe I should get. (i 
>> get something close but not 1.0. eg 0.96..)
>>
>> I am now wondering whether I am doubling up on the local_rhs in the 
>> assembling of cells on the boundary. To clarify, in step-20, f = 0 so that 
>> in the rhs_values, there is nothing added to the cells on the boundary 
>> other than that which is put in the boundary condition. for my problem, f 
>> is not equal to zero, so I am wondering whether I have extra contributions 
>> to the cells on the boundary and that's where the difference comes from. 
>>
>> Thank you for your help.
>>
>>
>>
>>
>>
>>
>>
>>
>> On Friday, February 22, 2019 at 10:51:11 PM UTC, Wolfgang Bangerth wrote:
>>>
>>> On 2/22/19 4:32 AM, jane...@jandj-ltd.com wrote: 
>>> > 
>>> > Wolfgang, so that's exactly what I had done with MMS, and that was 
>>> > verified, so i assumed the way i was imposing it was correct. 
>>>
>>> When you say "correct", you mean that you have two uncoupled equations? 
>>>
>>>
>>> > It then 
>>> > fails to get the correct value at the top. Initially, it's correct up 
>>> to 
>>> > about 4 decimal places, then to 2, then to the unit value. 
>>>
>>> Try to be precise. When you say "It then fails...", what has changed? 
>>> Above you say that things are correct, and now they are no longer. What 
>>> changes have you made in between? And what is "then to 4,... then to 2" 
>>> -- you are suggesting that something changes form one step to the next, 
>>> but what is it? 
>>>
>>> Best 
>>>   W. 
>>>
>>> -- 
>>>  
>>> Wolfgang Bangerth  email: bang...@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] Tutorial step-35 fails with my inp

[Wolfgang Bangerth]

Evgeny,

>    I found one problem with inp file - in line "1 0 quad    1 2 3 4". I 
> guess the second number should be nonzero here. I have try to set 1 and 
> get better results. Now I can see velocity changes around the 
> input-output bounds.

Is this an input file you created yourself, or are you saying that there 
is a mistake in the input files we distribute as part of step-35?


>    I think with cell face iterator I have to use cell quad iterator to 
> set all bounds when I create my inp, am I right?

I don't think I understand what you are trying to say here. An iterator 
can only point to a cell *or* a face. It can't be a "cell face iterator".


>    For input bound I use the same function as in the 35 tutorial. I 
> tried to set ConstantFunction() with value about 1 or -1 for the output 
> bound. But I stiil got the same error on the 3th step. I can see the 
> pressure was about 10^18 before that.

That suggests that you use boundary conditions that are not compatible 
with the incompressibility assumption. For example, if you prescribe 
velocity boundary values that point inward at every part of the 
boundary, you say that you want to have inflow everywhere -- but because 
the model assumes that your medium is incompressible and that there are 
no interior sinks for the fluids, this clearly can't work.


>    I get a good looking flow between input and output bounds when I use 
> ConstantFunction() with the same value for both input and output bounds. 
> But only for the about 5 steps. And then I can see the pressure grows 
> dramatically.
>    How I can chooze correct boundary conditions when I have variable 
> input boundary?

For one, you have to make sure that the boundary conditions are 
compatible with the model. In your case, this means that if you use 
Dirichlet or no-flux boundary conditions everywhere, then you need to 
have that
   \int_{partial Omega} n.u  = 0

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] ILU preconditioner and mesh-scale dependency?

[Wolfgang Bangerth]

On 3/13/19 8:18 PM, Jaekwang Kim wrote:
> 
> 
> I wonder now ...
> 
> 1. why smaller scale mesh only results such error,

You have two terms in your matrix, one from the advection term and one 
from the reaction term. If you scale the mesh by a factor of c, the 
height of the shape functions used in the reaction term does not change, 
but the gradient of shape functions used in the advection term is scaled 
by 1/c. As a consequence, the balance of the two terms that together 
make up the matrix changes and you apparently get into a situation where 
the resulting matrix is ill-conditioned or indeed not invertible.


> 2. If possible, is there any choice of better preconditionner for my 
> system other than SparseILU to go around this problem?

To make sure that the system you have assembled is correct, I usually 
first check with the SparseDirectUMFPACK solver. Once I know that the 
system is correct, I think about iterative solvers and preconditioners. 
In your case, the system you have is advection dominated, for which it 
is notoriously difficult to construct good preconditioners.

Best
  WB

-- 

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] Accuracy of Dirichlet condition for p in step-20

[Wolfgang Bangerth]

On 3/14/19 9:13 AM, jane@jandj-ltd.com wrote:
> Note that I think this was the issue.
> 
> You would need to put local_rhs = 0 within the part below to apply the 
> Dirichlet condition. Please do let me know if you think this is wrong 
> but I am getting the correct output. Thanks

That seems reasonable and correct.

Are you suggesting that step-20 could be made clearer by adding a line 
that in that program doesn't make a difference but would if it were used 
in a more general context? If so, could you suggest what to put and where?

Best & glad to hear that you solved the problem!
  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] what is the different between VectorTools::point_value and fe_values.get_function_values()?

[Wolfgang Bangerth]

On 3/14/19 1:28 AM, Phạm Ngọc Kiên wrote:
> 
> So for the first function 
> VectorTools::point_value(dof_handler,solution,point,values);
> I can specify the observed point in real coordiante, solution vector and 
> dof_handler to get the values, right?

Correct.


> For the second one, assuming that I am using QGauss for solving the 
> system, does it mean that if the point is on the edge of the cell then I 
> cannot get the solution at that point?

Yes. It can be any quadrature formula defined on the reference cell. It 
doesn't matter where the points are located.

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] Accuracy of Dirichlet condition for p in step-20

[jane . lee]

Note that I think this was the issue. 

You would need to put local_rhs = 0 within the part below to apply the 
Dirichlet condition. Please do let me know if you think this is wrong but I 
am getting the correct output. Thanks

for (unsigned int face_no=0;
face_no::faces_per_cell;
++face_no)
if (cell->face(face_no)->at_boundary()
&&
(cell->face(face_no)->boundary_id() == 1)) // top has boundary id 1 
On Thursday, March 14, 2019 at 2:24:55 PM UTC, jane...@jandj-ltd.com wrote:
>
> Hi Wolfgang, Sorry for the reply - for some reason i didn't get the 
> notification. 
>
> What i mean by correct is that I do have two uncoupled equations, but that 
> it is correct in that I had already verified it using what you suggested 
> (by using exact solutions - MMS).
>
> By it fails, I mean that the top pf value output is not equal to pr as it 
> should be, which I set in the weak form as in my first message in my 
> thread: 
>
> local_rhs(i) += -( fe_face_values_pf[velocities].value (i, q) *
> fe_face_values_pf.normal_vector(q) *
> pr_boundary_values[q] *
> fe_face_values_pf.JxW(q));
>
>  this is in the same way as step-20.
>
> I have now done further tests by just letting pr_boundary_values[q] to 1.0 
> in the blurb copied here - I do not get 1.0 as my value for the solution 
> for pf  on boundary with id 1, which I believe I should get. (i get 
> something close but not 1.0. eg 0.96..)
>
> I am now wondering whether I am doubling up on the local_rhs in the 
> assembling of cells on the boundary. To clarify, in step-20, f = 0 so that 
> in the rhs_values, there is nothing added to the cells on the boundary 
> other than that which is put in the boundary condition. for my problem, f 
> is not equal to zero, so I am wondering whether I have extra contributions 
> to the cells on the boundary and that's where the difference comes from. 
>
> Thank you for your help.
>
>
>
>
>
>
>
>
> On Friday, February 22, 2019 at 10:51:11 PM UTC, Wolfgang Bangerth wrote:
>>
>> On 2/22/19 4:32 AM, jane...@jandj-ltd.com wrote: 
>> > 
>> > Wolfgang, so that's exactly what I had done with MMS, and that was 
>> > verified, so i assumed the way i was imposing it was correct. 
>>
>> When you say "correct", you mean that you have two uncoupled equations? 
>>
>>
>> > It then 
>> > fails to get the correct value at the top. Initially, it's correct up 
>> to 
>> > about 4 decimal places, then to 2, then to the unit value. 
>>
>> Try to be precise. When you say "It then fails...", what has changed? 
>> Above you say that things are correct, and now they are no longer. What 
>> changes have you made in between? And what is "then to 4,... then to 2" 
>> -- you are suggesting that something changes form one step to the next, 
>> but what is it? 
>>
>> Best 
>>   W. 
>>
>> -- 
>>  
>> Wolfgang Bangerth  email: bang...@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] Accuracy of Dirichlet condition for p in step-20

[jane . lee]

Hi Wolfgang, Sorry for the reply - for some reason i didn't get the 
notification. 

What i mean by correct is that I do have two uncoupled equations, but that 
it is correct in that I had already verified it using what you suggested 
(by using exact solutions - MMS).

By it fails, I mean that the top pf value output is not equal to pr as it 
should be, which I set in the weak form as in my first message in my 
thread: 

local_rhs(i) += -( fe_face_values_pf[velocities].value (i, q) *
fe_face_values_pf.normal_vector(q) *
pr_boundary_values[q] *
fe_face_values_pf.JxW(q));

 this is in the same way as step-20.

I have now done further tests by just letting pr_boundary_values[q] to 1.0 
in the blurb copied here - I do not get 1.0 as my value for the solution 
for pf  on boundary with id 1, which I believe I should get. (i get 
something close but not 1.0. eg 0.96..)

I am now wondering whether I am doubling up on the local_rhs in the 
assembling of cells on the boundary. To clarify, in step-20, f = 0 so that 
in the rhs_values, there is nothing added to the cells on the boundary 
other than that which is put in the boundary condition. for my problem, f 
is not equal to zero, so I am wondering whether I have extra contributions 
to the cells on the boundary and that's where the difference comes from. 

Thank you for your help.








On Friday, February 22, 2019 at 10:51:11 PM UTC, Wolfgang Bangerth wrote:
>
> On 2/22/19 4:32 AM, jane...@jandj-ltd.com  wrote: 
> > 
> > Wolfgang, so that's exactly what I had done with MMS, and that was 
> > verified, so i assumed the way i was imposing it was correct. 
>
> When you say "correct", you mean that you have two uncoupled equations? 
>
>
> > It then 
> > fails to get the correct value at the top. Initially, it's correct up to 
> > about 4 decimal places, then to 2, then to the unit value. 
>
> Try to be precise. When you say "It then fails...", what has changed? 
> Above you say that things are correct, and now they are no longer. What 
> changes have you made in between? And what is "then to 4,... then to 2" 
> -- you are suggesting that something changes form one step to the next, 
> but what is it? 
>
> Best 
>   W. 
>
> -- 
>  
> Wolfgang Bangerth  email: bang...@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] what is the different between VectorTools::point_value and fe_values.get_function_values()?

[Jean-Paul Pelteret]

Dear Pham,

When you use the FEValuesExtractors in conjunction with the FEValues class 
using the operator[] like you’ve done here

> fe_values[vector_re]


then what is returned is an FEValuesViews type object. Since your vector_re 
extractor is of type FEValuesExtractors::vector, the returned object is a 
FEValuesViews::Vector. So the next thing to look at is the documentation for 
FEValuesViews::Vector::get_function_values() 
.
 It says that what is to be entered into the temp_solution_re vector  in this 
call
> fe_values[vector_re].get_function_values(solution, temp_solution_re);
is the vector-valued solution at all of the quadrature points on that cell (or, 
associated with the fe_values object). So, actually, the way that you 
initialise you temp_e_field_re vector should be
> std::vector> temp_e_field_re(quadrature.size());

and I wouldn’t be surprised if this code threw an error when running in debug 
mode, specifically when using QGauss of order > 2.

I hope that this helps you understand this a bit better.
Best,
Jean-Paul

> On 14 Mar 2019, at 08:28, Phạm Ngọc Kiên  wrote:
> 
> Thank you very much for your answer.
> 
> So for the first function  
> VectorTools::point_value(dof_handler,solution,point,values); 
> I can specify the observed point in real coordiante, solution vector and 
> dof_handler to get the values, right?
> 
> For the second one, assuming that I am using QGauss for solving the system, 
> does it mean that if the point is on the edge of the cell then I cannot get 
> the solution at that point?
> I am using the codes like this:
> const std::pair::active_cell_iterator, Point>
> cell_point = GridTools::find_active_cell_around_point(mapping, 
> dof_handler, point);
> 
> const double delta = 1e-20;
> Assert(GeometryInfo::distance_to_unit_cell(cell_point.second) < delta,
>ExcInternalError());
> 
> const Quadrature 
> quadrature(GeometryInfo::project_to_unit_cell(cell_point.second));
> FEValues fe_values(mapping, dof_handler.get_fe(), quadrature,
> update_values | update_gradients | 
> update_quadrature_points);
> fe_values.reinit(cell_point.first);
> 
> const FEValuesExtractors::Vector vector_re(0);
> const FEValuesExtractors::Vector vector_im(dim);
> 
> std::vector> temp_e_field_re(1);
> std::vector> temp_e_field_im(1);
> 
> //Get e_field at point
> fe_values[vector_re].get_function_values(solution, temp_solution_re);
> fe_values[vector_im].get_function_values(solution, temp_solution_im);
> 
> With this function, I see that it returns 9 components: 3 components for each 
> parameters: temp_solution_re[0], temp_solution_re[1], temp_solution_re[2] 
> when I print them out by:
> std::cout<<"result with get_function_values real: \t"<< temp_solution_re[0] 
> << ", "< std::cout<<"result with get_function_values imag: \t"<< temp_solution_im[0] 
> << ", "< 
> I really want to know which one is the solution for my vector in x,y,z 
> direction?
> 
> Best regards,
> Pham Ngoc Kien
> 
> Vào 14:08:11 UTC+9 Thứ Năm, ngày 14 tháng 3 năm 2019, Wolfgang Bangerth đã 
> viết:
> On 3/13/19 8:21 PM, Phạm Ngọc Kiên wrote: 
> > I am testing my codes for output the solution at a point in my numerical 
> > model. 
> > I saw in the library 2 ways to do this task. The first one is to use 
> > VectorTools::point_value() function and the second one is 
> > fe_values.get_function_values(). 
> 
> These functions are really used in very different contexts. The first of 
> these 
> is when you want to evaluate a solution vector at some arbitrary point in the 
> domain. As a consequence, you have to specify this point when you call the 
> function. 
> 
> The second function is used when you want to evaluate the function at very 
> specific points, namely the quadrature points of the current cell. You can of 
> course only use this if the point at which you want to evaluate the solution 
> happens to be a quadrature point of the current cell, and not just any point 
> anywhere. In return the second function is then vastly more efficient than 
> the 
> first. 
> 
> Does this help? 
> 
> Best 
>   W. 
> 
> -- 
>  
> Wolfgang Bangerth  email: bang...@ <>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...@googlegro

Re: [deal.II] what is the different between VectorTools::point_value and fe_values.get_function_values()?

[Phạm Ngọc Kiên]

Thank you very much for your answer.

So for the first function  
VectorTools::point_value(dof_handler,solution,point,values); 

I can specify the observed point in real coordiante, solution vector and 
dof_handler to get the values, right?

For the second one, assuming that I am using QGauss for solving the system, 
does it mean that if the point is on the edge of the cell then I cannot get 
the solution at that point?
I am using the codes like this:

const std::pair::active_cell_iterator, Point>
cell_point = GridTools::find_active_cell_around_point(mapping, 
dof_handler, point);

const double delta = 1e-20;
Assert(GeometryInfo::distance_to_unit_cell(cell_point.second) < delta,
   ExcInternalError());

const Quadrature 
quadrature(GeometryInfo::project_to_unit_cell(cell_point.second));
FEValues fe_values(mapping, dof_handler.get_fe(), quadrature,
update_values | update_gradients | 
update_quadrature_points);
fe_values.reinit(cell_point.first);

const FEValuesExtractors::Vector vector_re(0);
const FEValuesExtractors::Vector vector_im(dim);

std::vector> temp_e_field_re(1);
std::vector> temp_e_field_im(1);

//Get e_field at point
fe_values[vector_re].get_function_values(solution, temp_solution_re);
fe_values[vector_im].get_function_values(solution, temp_solution_im);


With this function, I see that it returns 9 components: 3 components for 
each parameters: temp_solution_re[0], temp_solution_re[1], temp_solution_re[
2] 
when I print them out by:

std::cout<<"result with get_function_values real: \t"<< temp_solution_re[0] << 
", "<
> On 3/13/19 8:21 PM, Phạm Ngọc Kiên wrote: 
> > I am testing my codes for output the solution at a point in my numerical 
> model. 
> > I saw in the library 2 ways to do this task. The first one is to use 
> > VectorTools::point_value() function and the second one is 
> > fe_values.get_function_values(). 
>
> These functions are really used in very different contexts. The first of 
> these 
> is when you want to evaluate a solution vector at some arbitrary point in 
> the 
> domain. As a consequence, you have to specify this point when you call the 
> function. 
>
> The second function is used when you want to evaluate the function at very 
> specific points, namely the quadrature points of the current cell. You can 
> of 
> course only use this if the point at which you want to evaluate the 
> solution 
> happens to be a quadrature point of the current cell, and not just any 
> point 
> anywhere. In return the second function is then vastly more efficient than 
> the 
> first. 
>
> Does this help? 
>
> Best 
>   W. 
>
> -- 
>  
> Wolfgang Bangerth  email: bang...@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] what is the different between VectorTools::point_value and fe_values.get_function_values()?

[Phạm Ngọc Kiên]

Thank you very much for your answer.
So for the first function
VectorTools::point_value(dof_handler,solution,point,values);

I can specify the observed point in real coordiante, solution vector and
dof_handler to get the values right?

For the second one, assuming that I am using QGauss for solving the system,
can I use QGaussLobatto to get the quadrature points on the edge of a cell
and get the values from this?
With the second function, I see that it returns 9 components: 3 for each
parameters: temp_solution_re[0], temp_solution_re[1], temp_solution_re[3]
when I call fe_values[vector_re].get_function_values(solution,
temp_solution_re);
I really want to know which one is the solution for my vector in x,y,z
direction?

Best regards,
Pham Ngoc Kien

Vào Th 5, 14 thg 3, 2019 vào lúc 14:08 Wolfgang Bangerth <
bange...@colostate.edu> đã viết:

> On 3/13/19 8:21 PM, Phạm Ngọc Kiên wrote:
> > I am testing my codes for output the solution at a point in my numerical
> model.
> > I saw in the library 2 ways to do this task. The first one is to use
> > VectorTools::point_value() function and the second one is
> > fe_values.get_function_values().
>
> These functions are really used in very different contexts. The first of
> these
> is when you want to evaluate a solution vector at some arbitrary point in
> the
> domain. As a consequence, you have to specify this point when you call the
> function.
>
> The second function is used when you want to evaluate the function at very
> specific points, namely the quadrature points of the current cell. You can
> of
> course only use this if the point at which you want to evaluate the
> solution
> happens to be a quadrature point of the current cell, and not just any
> point
> anywhere. In return the second function is then vastly more efficient than
> the
> first.
>
> Does this help?
>
> 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.