Hi,
sorry for asking your help again, that I really appreciated.
Today's issue is about imposing boundary conditions for a three-fields
formulation.
In brief, I defined a 4 nodes square element, of type:
fe(FE_Q<dim>(poly_degree), dim, // displacement
FE_DGPMonomial<dim>(poly_degree - 1), 1, // pressure
FE_DGPMonomial<dim>(poly_degree - 1), 1), // dilatation
with poly_degree = 1.
While debugging on a single element, which has 10 dofs in 2D ( 4x2
displacements + 1 pressure + 1 dilatation), I aim at imposing Dirichlet
boundary conditions on displacements, whereas body forces are zero.
The cell matrix and the cell_rhs are the following.
*Cell matrix = *
*{ 0.7037, 0.0278, -0.2037, -0.4722, -0.1481, 0.4722, -0.3519,
-0.0278, -1.0000, 0.0000}, *
*{ 0.0278, 0.7037, 0.4722, -0.1481, -0.4722, -0.2037, -0.0278,
-0.3519, -1.0000, 0.0000}, *
*{ -0.2037, 0.4722, 0.7037, -0.0278, -0.3519, 0.0278, -0.1481,
-0.4722, 1.0000, 0.0000}, *
*{ -0.4722, -0.1481, -0.0278, 0.7037, 0.0278, -0.3519, 0.4722,
-0.2037, -1.0000, 0.0000}, *
*{ -0.1481, -0.4722, -0.3519, 0.0278, 0.7037, -0.0278, -0.2037,
0.4722, -1.0000, 0.0000}, *
*{ 0.4722, -0.2037, 0.0278, -0.3519, -0.0278, 0.7037, -0.4722,
-0.1481, 1.0000, 0.0000}, *
*{ -0.3519, -0.0278, -0.1481, 0.4722, -0.2037, -0.4722, 0.7037,
0.0278, 1.0000, 0.0000}, *
*{ -0.0278, -0.3519, -0.4722, -0.2037, 0.4722, -0.1481, 0.0278,
0.7037, 1.0000, 0.0000}, *
*{ -1.0000, -1.0000, 1.0000, -1.0000, -1.0000, 1.0000, 1.0000,
1.0000, 0.0000, -4.0000}, *
*{ 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000,
0.0000, -4.0000, 6.6667}, *
*Cell rhs =0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00
0.000e+00 0.000e+00 0.000e+00 0.000e+00 *
They seem to be correct. After calling
system_matrix.compress(VectorOperation::add);
system_rhs.compress(VectorOperation::add);
the system matrix coincides with the cell matrix, the same for the rhs. To
impose BCs I followed the scheme of step-18, as follows:
std::map<types::global_dof_index,double> boundary_values;
ComponentMask displacements_mask = fe.component_mask (displacements);
VectorTools::interpolate_boundary_values (dof_handler,
0,
IncrementalDirichletBoundaryValues<dim>( TIDM, NRit, GPDofs ),
boundary_values,
displacements_mask);
PETScWrappers::MPI::Vector tmp (locally_owned_dofs,mpi_communicator);
MatrixTools::apply_boundary_values (boundary_values,
system_matrix,
tmp,
system_rhs,
false);
incremental_displacement = tmp;
where the class IncrementalDirichletBoundaryValues implements the Dirichlet
BC and GPDofs = 4.
The boundary values after VectorTools::interpolate_boundary_values turn out
to be:
*boundary_values = -0.0500, -0.0500, -0.0500, 0.0500, 0.0500,
-0.0500, 0.0500, 0.0500, 0.0000, 0.0000, *
*tmp = -0.0500, -0.0500, -0.0500, 0.0500, 0.0500, -0.0500, 0.0500,
0.0500, 0.0000, 0.0000, *
which are pretty much what I expected. However, after MatrixTools
::apply_boundary_values the system matrix and the rhs are as follows,
* System matrix = *
*{ 0.7037, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000,
0.0000, 0.0000, 0.0000}, *
*{ 0.0000, 0.7037, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000,
0.0000, 0.0000, 0.0000}, *
*{ 0.0000, 0.0000, 0.7037, 0.0000, 0.0000, 0.0000, 0.0000,
0.0000, 0.0000, 0.0000}, *
*{ 0.0000, 0.0000, 0.0000, 0.7037, 0.0000, 0.0000, 0.0000,
0.0000, 0.0000, 0.0000}, *
*{ 0.0000, 0.0000, 0.0000, 0.0000, 0.7037, 0.0000, 0.0000,
0.0000, 0.0000, 0.0000}, *
*{ 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.7037, 0.0000,
0.0000, 0.0000, 0.0000}, *
*{ 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.7037,
0.0000, 0.0000, 0.0000}, *
*{ 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000,
0.7037, 0.0000, 0.0000}, *
*{ -1.0000, -1.0000, 1.0000, -1.0000, -1.0000, 1.0000, 1.0000,
1.0000, 0.0000, -4.0000}, *
*{ 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000,
0.0000, -4.0000, 6.6667}, *
*System_rhs = -0.0352, -0.0352, -0.0352, 0.0352, 0.0352, -0.0352,
0.0352, 0.0352, 0.0000, 0.0000, *
The* System_rhs *seems correct, but the matrix, particularly the 9th row,
does not, but I am not sure why.
Many thanks in advance,
Alberto
--
Informativa sulla Privacy: http://www.unibs.it/node/8155
--
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.
