Issue with Applying Dirichlet Boundary Conditions

[Alex]

^{Hello everyone,}

^{I am modifying the Step-44 (Quasi-Static Finite-Strain Compressible Elasticity) tutorial in deal.II to incorporate thermo-mechanical equations. However, I am facing an issue with applying temperature Dirichlet boundary conditions.}

^{Problem Description:}

^{I am applying the temperature boundary condition using:}

^{VectorTools::interpolate_boundary_values(}

    ^{dof_handler_ref,}

    ^boundary_id,

    ^{Functions::ConstantFunction<dim>(current_applied_temperature_0, dim + 1),} 

    ^constraints,

    ^{fe.component_mask(Temperature)}

^);

^{After that, before solving the system I apply:}

^{constraints.condense(tangent_matrix, system_rhs);}

^{Despite this, my solution converges, but the computed nodal temperatures are negative and completely incorrect.}

^{However, when I apply boundary conditions using:}

^{MatrixTools::apply_boundary_values(boundary_values,tangent_matrix,solution_n,system_rhs);}

^{I obtain the correct nodal temperatures.}

 Could someone clarify what am I doing wrong in my first approach?  Thanks for your support

 

 

 

 

 

[Alex]

Hi all,

I would like to add a few statements to my above question, I have tried to boil down the issue and noticed that the issue arises with the global right-hand side (RHS) after applying AffineConstraints::condense(). Specifically, my RHS is initially zero, but after applying constraints, unexpected large negative values appear at certain degrees of freedom (DoFs).

I have also attached my code below.

Setup

• I have 2 elements with 6 nodes.
• The far-left (200 kelvin) and far-right nodes  (0 Kelvin) have prescribed temperature Dirichlet boundary conditions, while the middle nodes remain unconstrained.
• The element-wise RHS (cell_rhs) remains zero before assembling the global system.
• I print the global RHS before and after constraints for debugging.

Observed Behavior

At time step = 1, and Newton iteration = 0, I observe the following outputs:

1. Before Applying Constraints:

Global RHS BEFORE Constraints:

0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00 

2. After Applying Constraints:

Global RHS AFTER Constraints:

0.000e+00  0.000e+00  -1.900e+04  -1.900e+04  0.000e+00  0.000e+00 

Can someone guide me to where am I going wrong?

Thanks for your support.

Thanks for your support !

 

 

 

 

 

--
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.
To view this discussion visit https://groups.google.com/d/msgid/dealii/0816bcff-fef9-40c9-972a-4ca7f626e75bn%40googlegroups.com.

[Jean-Paul Pelteret]

Hi Alex,

I haven't inspected your code in detail, but I would suggest that the problem might be related to which system you are applying the constraints, and into which vector you are distributing the constraints. It looks to me that you are trying to specify the total displacement and temperature at some locations, 

      const double current_applied_Displacement = total_displacement *time_steps;
      const double current_applied_temperature_0 = total_temperature_0*time_steps;
      double current_applied_temperature_L= total_temperature_L*time_steps;

while the constraints are being distributed to the solution increment

constraints.distribute(newton_update);

You might want to try either reformulating your boundary conditions as updates to the BCs applied at the previous timestep (i.e. constrain the increment) or you should refactor the code (assembly, constraints creation + something before the Newton iterations start) such that you apply the constraints directly to the solution vector itself.

I hope that this helps.

Best,

Jean-Paul

To view this discussion visit https://groups.google.com/d/msgid/dealii/CAK1JoDy69EbXrjK-aYx5DjBCbxdaFFEjWdeyQy_PO79601sd4Q%40mail.gmail.com.

[Alex]

Dear Jean-Paul,

Thank you for your response—I really appreciate your insights.

Regarding your suggestion: "you should refactor the code (assembly, constraints creation + something before the Newton iterations start) such that you apply the constraints directly to the solution vector itself."

When I apply my temperature boundary conditions differently using:

MatrixTools::apply_boundary_values(boundary_values, tangent_matrix, solution_n, system_rhs);

(which modifies solution_n instead of newton_update) before starting the Newton-Raphson iterations, everything works as expected, and I achieve the correct results.

However, when I instead use:

VectorTools::interpolate_boundary_values(
    dof_handler_ref,
    boundary_id,
    Functions::ConstantFunction<dim>(current_applied_temperature_0, dim + 1),
    constraints,
    fe.component_mask(Temperature)
);

  copy_local_to_global_ASM(const PerTaskData_ASM &data)

  {

  const AffineConstraints<double> &constraints = data.solid->constraints;

  BlockSparseMatrix<double> &tangent_matrix = const_cast<Solid<dim,NumberType> *>(data.solid)->tangent_matrix;

  BlockVector<double> &system_rhs = const_cast<Solid<dim,NumberType> *>(data.solid)->system_rhs;

 

constraints.condense(tangent_matrix, system_rhs);

}

 constraints.distribute(newton_update);

I notice a discrepancy—though the Newton-Raphson iterations still converge, the solution does not fully align with expectations. This seems quite odd, given that both approaches should theoretically produce the same results.

Would you have any thoughts on why this might be happening? Any guidance would be greatly appreciated.

 Regards,

Pradeep.

[Wolfgang Bangerth]

On 3/19/25 00:51, Jean-Paul Pelteret wrote:
>
> You might want to try either reformulating your boundary conditions as updates
> to the BCs applied at the previous timestep (i.e. constrain the increment) or
> you should refactor the code (assembly, constraints creation + something
> before the Newton iterations start) such that you apply the constraints
> directly to the solution vector itself.

There is also a discussion on boundary values for Newton methods in step-15
that you might want to read through. In short, you likely only want to apply
Dirichlet boundary conditions to the first Newton iteration, and then use zero
boundary conditions after that.

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Alex]

Dear Prof. Wolfgang,

Thank you very much for your kind response.

I sincerely appreciate your insights. However, I believe this aspect has already been addressed in step-44 through the conditions mentioned below. 

    // Since the constraints are different at different Newton iterations, we

    // need to clear the constraints matrix and completely rebuild

    // it. However, after the first iteration, the constraints remain the same

    // and we can simply skip the rebuilding step if we do not clear it.

    if (it_nr > 1)

      return;

    const bool apply_dirichlet_bc = (it_nr == 0);

To verify this, I printed the values of the constraints at each iteration. Based on my observations, the constraints are populated with values only at Iteration "0," and in all subsequent iterations, they remain zero.

Please let me know if you have any thoughts or if I might have overlooked anything. I truly appreciate your guidance.

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

To view this discussion visit https://groups.google.com/d/msgid/dealii/d229e447-a652-4877-8b03-5659c14b83f5%40colostate.edu.

[Alex]

Hi Everyone,

I hope you're all doing well.

I was wondering if anyone might have some insight into why this isn't working. I've also reviewed the videos on applying non-Dirichlet boundary conditions in deal.II, and everything seems to align with the code.

Despite that, the issue still persists!

Thank you in advance for your time and help.

[Wolfgang Bangerth]

Alex,

This is the starting point on how to debug this. If I understand you
right, you expect the constraints object to contain things also at
iterations other than the first one, but it does not. You you have
something that seems wrong to you. Now work backward to find out why the
object does not contain anything. Where is the place where it's supposed
to be filled? Why is that not happening?

Best
W.