Dear all,
I am a master student working working with Glacial Isostatic Adjustment at
Uppsala University, and I find the deal.ii a very exciting and powerful
tool.
Basically I have arranged for my code the stokes flow (step-22) for the
model in use. The idea is very simple and I believe it can be achieved by
following the ideas in step-22.
I would like to ask a question that bugs me, the block structure I use
(similar as in step-22) is of the form:
| A Bt |
| B -C |
where Bt is the transpose of B and C is the mass matrix for the pressure.
The issue is that, even though I use
DoFRenumbering::component_wise (dof_handler, block_component);
the final global matrix gets the desired block structure, but the local
matrix does not. My first question is: *how to make the local matrix have
this block pattern as well?*
I need to reorder the local matrix so that the indexes are of the form:
( u0, u1, u2,...u8, v0, v1, v2, ... v8, p0, p1, p2, p3 )
The way the local indexes are in step-22 is:
(u0, v0, p0, u1, v1, p1, ... ,u3 ,v3 ,p3 ,u4 ,v4 ,u5 ,v5 ,u6 ,v6 ... u8
,v8)
where u = horizontal velocity, v = vertical velocity, p = pressure.
The transformation I do is rather artificial and it does not look good
(prone to errors):
...
const unsigned int dim_disp = 9, // Q2 nodes displacements
dim_u = dim*dim_disp, // dim * Q2 nodes
dim_p = 4; // 1 * Q1 nodes
unsigned int l_p[] = {2,5,8,11};
unsigned int order[] = {0,3,6,9,12,14,16,18,20, 1,4,7,10,13,15,17,19,21,
2,5,8,11};
FullMatrix<double> l_A (dim_u,dim_u), // local matrices for the
element-by-element Schur construction, l_A stands for local
l_Bt (dim_u,dim_p),
l_B (dim_p,dim_u),
l_C (dim_p,dim_p),
l_S (dim_p,dim_p),
aux (dim_u,dim_u);
...
and then I transform the cell_matrix obtained by the regular way:
...
// extract here using velocities, pressure reordering
for (unsigned int i=0; i<dofs_per_cell; ++i){ // shape index i
for (unsigned int j=0; j < dofs_per_cell; ++j){ // shape index j
cell_ordered(i,j) = cell_matrix(order[i],order[j]); // local
matrix ordered by u0,u1...un,v0,v1...vn,p0,p1...pm
if( i < dim_u ){
if(j < dim_u){
l_A(i,j) = cell_ordered(i,j); // i < dim_u, j < dim_u
aux(i,j) = l_A(i,j);
}else{
l_Bt(i,j-dim_u) = cell_ordered(i,j); // i < dim_u, j >=
dim_u
}
}else{
if(j < dim_u){
l_B(i-dim_u,j) = cell_ordered(i,j); // i >= dim_u, j <
dim_u
}else{
l_C(i-dim_u,j-dim_u) = cell_ordered(i,j); // i >=
dim_u, j >= dim_u
}
}
}
}
...
After building the local Schur (l_S), I need to reasemble in the old index
form (instead of C) to be able to assemble local-to-global:
...
for (unsigned int i=0; i< dim_p; ++i){ // shape index i
for (unsigned int j=0; j < dim_p; ++j){ // shape index j
cell_precond(l_p[i],l_p[j]) = l_S(i,j);
}// j
}// i
...
and then distribute local to global ....
...
// Local to global
cell->get_dof_indices (local_dof_indices);
constraints.distribute_local_to_global (cell_matrix, cell_rhs,
local_dof_indices,
system_matrix, system_rhs);
constraints.distribute_local_to_global (cell_precond,
local_dof_indices,
system_preconditioner);
...
I believe it is too much hassle and if there is no way to achieve this it
would be nice to include it.
another question that I have is *how to get the grid size of the
cell?* something
like dx, or the diagonal h = sqrt(dx^2+dy^2+dz^2).
I hope to hear from you,
--
Kind regards,
Juan Carlos Araújo-Cabarcas.
Uppsala University
_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
