Getting RHS values at nodes with DBC

[RAJAT ARORA]

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.

[Jean-Paul Pelteret]

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<dof_handler.n_dofs(); ++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

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

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

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

[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: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[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

[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 :-)

[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:

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?

[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

[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