Re: [deal.II] Getting RHS values at nodes with DBC

[Giorgos Kourakos]

As always there is a "deal.ii" way of doing the calculations.

The FEValues::get_function_gradients() is actually what I needed.

I got a little bit confused with the FEValues::get_normal_vectors(). This 
seems to return the normal direction of the face and doesn't give any 
further info about the flow direction.

Thank you for your help.

Giorgos



On Tuesday, September 3, 2019 at 8:37:01 AM UTC-7, Daniel Arndt wrote:
>
> Giorgos,
>
> You can find a discussion for how constraints are handled in deal.II in 
> https://www.dealii.org/developer/doxygen/deal.II/group__constraints.html.
>
> Thank you for your feedback.
>> I'll try to explain what I try to calculate with the help of a standard 
>> FEM textbook:
>>
>> [...]
>>
>> In the book the global system is partitioned so that the first row 
>> contains the known U1 dofs and the second row the unknown the dofs U2.
>> They solve the system for the U2 and then they calculate F1.
>>
>
> In short: we normally don't condense the linear system but use a diagonal 
> matrix for the constrained rows instead (and modify the right-hand side 
> accordingly).
>  
>
>> I think that the F1 vector  contains, among others, the fluxes from the 
>> imposed BCs.
>>
> In case you only have Dirichlet boundary conditions, K11 is diagonal, K12 
> is zero and F1 contains the (possibly scaled) Dirichlet boundary values.
>
> Could you clarify in mathematical terms what you mean by "fluxes"? Maybe, 
> we are just using different terms fo the same thing.
>
> In case you are interested in n\cdot In u at the boundary, you can use 
> FEValues::get_function_gradients() and FEValues::get_normal_vectors() with 
> a quadrature formel that contains
> the points you are interested in on a cell in reference coordinates. This 
> is the approach VectorTools::point_gradient() is using for returning the 
> gradient for a single point in real coordinates.
>
> 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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/a7a0cfa9-1a31-4133-bad0-9e014d88b4d5%40googlegroups.com.

Re: [deal.II] Getting RHS values at nodes with DBC

[Daniel Arndt]

Giorgos,

You can find a discussion for how constraints are handled in deal.II in
https://www.dealii.org/developer/doxygen/deal.II/group__constraints.html.

Thank you for your feedback.
> I'll try to explain what I try to calculate with the help of a standard
> FEM textbook:
>
> [...]
>
> In the book the global system is partitioned so that the first row
> contains the known U1 dofs and the second row the unknown the dofs U2.
> They solve the system for the U2 and then they calculate F1.
>

In short: we normally don't condense the linear system but use a diagonal
matrix for the constrained rows instead (and modify the right-hand side
accordingly).


> I think that the F1 vector  contains, among others, the fluxes from the
> imposed BCs.
>
In case you only have Dirichlet boundary conditions, K11 is diagonal, K12
is zero and F1 contains the (possibly scaled) Dirichlet boundary values.

Could you clarify in mathematical terms what you mean by "fluxes"? Maybe,
we are just using different terms fo the same thing.

In case you are interested in n\cdot In u at the boundary, you can use
FEValues::get_function_gradients() and FEValues::get_normal_vectors() with
a quadrature formel that contains
the points you are interested in on a cell in reference coordinates. This
is the approach VectorTools::point_gradient() is using for returning the
gradient for a single point in real coordinates.

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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/CAOYDWb%2B0Vj_LTboMTc2Rqo0%2BMLc7VQpNK0QfZFgeqabk_f7qbg%40mail.gmail.com.

Re: [deal.II] Getting RHS values at nodes with DBC

[Giorgos Kourakos]

Thank you for your feedback.
I'll try to explain what I try to calculate with the help of a standard FEM 
textbook:

[image: CaptureFEMbook.PNG]

In the book the global system is partitioned so that the first row contains 
the known U1 dofs and the second row the unknown the dofs U2.
They solve the system for the U2 and then they calculate F1.

I think that the F1 vector  contains, among others, the fluxes from the 
imposed BCs.

Does this make sense?

Can I achieve something like that with deal.ii?

-- 
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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/5fda3e26-e4bc-476a-9929-728a44858a20%40googlegroups.com.

Re: [deal.II] Getting RHS values at nodes with DBC

[Wolfgang Bangerth]

> In my case I solve a simple laplace problem, that represent groundwater flow 
> with non zero dirichlet boundary conditions
> My goal is to estimate how much water comes in and out from the dirichlet BC 
> nodes.

You're asking the wrong question. What you are looking for is a *flux*, which 
is an integral over part of the boundary. That's how you should compute it -- 
which implies you'll want to use FEFaceValues as always when face integrals 
are involved. Specifically, the flux through a single point is of course zero.


> The way I could do this, is after I solve the reduced version of the full 
> linear system A*U = F
> to calculate
> DBC = A*U-F

I'm not clear what you're trying to do here, but if A*U=F, then DBC is of 
course zero. I suspect that you mean for either A, U, or F to be something 
slightly different than what you had in the linear system above, but it's not 
clear to me -- though it may also not matter in view of what I mentioned above 
:-)

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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/1abaa9aa-1464-28e0-22fc-56d5e379fa4e%40colostate.edu.

Re: [deal.II] Getting RHS values at nodes with DBC

[Giorgos Kourakos]

This is a year old post but I'm dealing with the same issue so I decided to 
continue here instead creating a new one.

In my case I solve a simple laplace problem, that represent groundwater 
flow with non zero dirichlet boundary conditions
My goal is to estimate how much water comes in and out from the dirichlet 
BC nodes.

The way I could do this, is after I solve the reduced version of the full 
linear system A*U = F 
to calculate
DBC = A*U-F

When I call 
system_rhs.reinit (locally_owned_dofs, mpi_communicator);
system_matrix.reinit (locally_owned_dofs,
  locally_owned_dofs,
  dynamic_sparsity_pattern,
  mpi_communicator);


to initialize the A and F, are these the reduced (without the dofs of DBC) 
or the full ones?
  
is
system_matrix.vmult(output, system_rhs);

equivalent to ?
A*U (where A includes the rows and columns of dirichlet BC)

Thank you

-- 
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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/9d205351-23b1-4dbb-aebf-581ea3b19f63%40googlegroups.com.

Re: [deal.II] Getting RHS values at nodes with DBC

[Wolfgang Bangerth]

On 08/01/2018 03:34 PM, RAJAT ARORA wrote:


If I understand correctly, I dont think that will work. The reason is 
because of the Dirichlet boundary conditions, the value of the rhs at those 
nodes is unknown until you solve for the system.


I think it's unclear to me what exactly you want in these entries. What you 
suggest here is that in entry 'i' of a rhs vector (where 'i' is a boundary 
DoF) you want the correct value U_i of the solution at this node.


But in general, the 'i'th entry of the rhs vector (for any other 'i') contains 
the integral

  rhs_i   =  \int  \phi_i(x)  f(x)  dx

So can you describe, in mathematical terms, what exactly you are looking for? 
Or, maybe even better, what you want to do with this information?


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] Getting RHS values at nodes with DBC

[Jean-Paul Pelteret]

> However, now how do I get the value of the rhs at the nodes where DBC were 
> specified?


So, are you wanting the values to which the DoFs on the Dirichlet boundary are 
going to be prescribed? This is the same as the values stored in the constraint 
matrix (assuming that it holds only Dirichlet constraints). You could use 
DoFTools::extract_boundary_dofs() to find the DoF indices of those that are 
constrained (i.e. on the specified Dirichlet boundaries), and then extract the 
inhomogeneity associated with that line in the constraint matrix.

> On 01 Aug 2018, at 23:34, RAJAT ARORA  wrote:
> 
> Hi Jean,
> 
> Thanks for your reply. 
> 
> If I understand correctly, I dont think that will work. The reason is because 
> of the Dirichlet boundary conditions, the value of the rhs at those nodes is 
> unknown until you solve for the system.
> 
> 
> 
> 
> On Wednesday, August 1, 2018 at 10:08:12 AM UTC-4, Jean-Paul Pelteret wrote:
> Hi Rajat,
> 
> Could you perhaps assemble a second RHS vector with an empty constraint 
> matrix, and then extract the entries in this vector that are constrained by 
> the primary constraint matrix? i.e. something like
> 
> PETScVector unconstrained_rhs_vector;
> ConstraintMatrix empty_constraints;
> empty_constraints.close();
> 
> // Assembly
> empty_constraints.distribute_local_to_global(local_rhs_vector, 
> local_dof_indices, unconstrained_rhs_vector);
> 
> …..
> // Extraction
> 
> for (unsigned int d=0; d   if (locally_owned_dofs.is_element(d) && constraints.is 
> _constrained(d) && 
> !hanging_node_constraints.is_constrained(d))
> const double value = unconstrained_rhs_vector(d);
> 
> Best,
> Jean-Paul
> 
>> On 01 Aug 2018, at 05:33, RAJAT ARORA gmail.com 
>> > wrote:
>> 
>> Hi,
>> 
>> I am using deal.ii linked with MPI and Petsc in my PhD work. 
>> 
>> I have a problem to which I have a temporary solution but I think it can be 
>> done in a better way.
>> 
>> The problem is as follows: Lets say we are solving an equation with some 
>> Dirichlet boundary conditions (DBC).
>> As I assemble the system (using constraints.distribute_local_to_global), the 
>> matrix  and rhs are modified to give zero values for the solution at the 
>> dofs where DBC are specified.
>> 
>> So, I call constraints.distribute to correctly set the known values.
>> 
>> However, now how do I get the value of the rhs at the nodes where DBC were 
>> specified?
>> 
>> The way I do this is to simultaneously assemble another stiffness matrix 
>> without any knowledge of Dirichlet constraints. I then multiply the obtained 
>> solution with this stiffness matrix
>> to the full known rhs.
>> 
>> Is this optimal way of doing this or can this be done in a better way 
>> (without assembly the matrix twice)? 
>> 
>> The other way I can think of doing the same thing but on a local element 
>> level. If I donot need it everywhere then this may be cheaper.
>> 
>> 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+un...@ <>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 
> .

-- 
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] Getting RHS values at nodes with DBC

[RAJAT ARORA]

Hi Jean,

Thanks for your reply. 

If I understand correctly, I dont think that will work. The reason is 
because of the Dirichlet boundary conditions, the value of the rhs at those 
nodes is unknown until you solve for the system.




On Wednesday, August 1, 2018 at 10:08:12 AM UTC-4, Jean-Paul Pelteret wrote:
>
> Hi Rajat,
>
> Could you perhaps assemble a second RHS vector with an empty constraint 
> matrix, and then extract the entries in this vector that are constrained by 
> the primary constraint matrix? i.e. something like
>
> PETScVector unconstrained_rhs_vector;
> ConstraintMatrix empty_constraints;
> empty_constraints.close();
>
> // Assembly
> empty_constraints.distribute_local_to_global(local_rhs_vector, 
> local_dof_indices, unconstrained_rhs_vector);
>
> …..
> // Extraction
>
> for (unsigned int d=0; d   if (locally_owned_dofs.is_element(d) && constraints.is_constrained(d) 
> && !hanging_node_constraints.is_constrained(d))
> const double value = unconstrained_rhs_vector(d);
>
> Best,
> Jean-Paul
>
> On 01 Aug 2018, at 05:33, RAJAT ARORA > 
> wrote:
>
> Hi,
>
> I am using deal.ii linked with MPI and Petsc in my PhD work. 
>
> I have a problem to which I have a temporary solution but I think it can 
> be done in a better way.
>
> The problem is as follows: Lets say we are solving an equation with some 
> Dirichlet boundary conditions (DBC).
> As I assemble the system (using constraints.distribute_local_to_global), 
> the matrix  and rhs are modified to give zero values for the solution at 
> the dofs where DBC are specified.
>
> So, I call constraints.distribute to correctly set the known values.
>
> However, now how do I get the value of the rhs at the nodes where DBC were 
> specified?
>
> The way I do this is to simultaneously assemble another stiffness matrix 
> without any knowledge of Dirichlet constraints. I then multiply the 
> obtained solution with this stiffness matrix
> to the full known rhs.
>
> Is this optimal way of doing this or can this be done in a better way 
> (without assembly the matrix twice)? 
>
> The other way I can think of doing the same thing but on a local element 
> level. If I donot need it everywhere then this may be cheaper.
>
> 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+un...@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.