From: kanschat Date: Wed, 14 Mar 2012 17:30:50 +0000 (+0000) Subject: local blocks for Assembler::MGMatrixSimple X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=e327d26b963937a26fb51e89f011a2c971a9de27;p=dealii-svn.git local blocks for Assembler::MGMatrixSimple git-svn-id: https://svn.dealii.org/trunk@25280 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/include/deal.II/meshworker/dof_info.h b/deal.II/include/deal.II/meshworker/dof_info.h index 40addebb6c..f09b6ab2d3 100644 --- a/deal.II/include/deal.II/meshworker/dof_info.h +++ b/deal.II/include/deal.II/meshworker/dof_info.h @@ -187,7 +187,9 @@ namespace MeshWorker * DoFInfoBox. */ DoFInfo (); - + + /// Set up local block indices + void set_block_indices (); /// Fill index vector with active indices void get_indices(const typename DoFHandler::cell_iterator& c); diff --git a/deal.II/include/deal.II/meshworker/dof_info.templates.h b/deal.II/include/deal.II/meshworker/dof_info.templates.h index 76faf577b6..3ca29e4bc2 100644 --- a/deal.II/include/deal.II/meshworker/dof_info.templates.h +++ b/deal.II/include/deal.II/meshworker/dof_info.templates.h @@ -33,7 +33,24 @@ namespace MeshWorker DoFInfo::DoFInfo() {} + + template + void + DoFInfo::set_block_indices() + { + for (unsigned int i=0;i + bi = block_info->local().global_to_local(this->block_info->renumber(i)); + indices_by_block[bi.first][bi.second] = indices_org[i]; + } + // Remove this after + // changing block codes + for (unsigned int i=0;iblock_info->renumber(i)] = indices_org[i]; + } + template void DoFInfo::get_indices(const typename DoFHandler::cell_iterator& c) @@ -48,16 +65,7 @@ namespace MeshWorker { indices_org.resize(c->get_fe().dofs_per_cell); c->get_dof_indices(indices_org); - for (unsigned int i=0;i - bi = block_info->local().global_to_local(this->block_info->renumber(i)); - indices_by_block[bi.first][bi.second] = indices_org[i]; - } - // Remove this after - // changing block codes - for (unsigned int i=0;iblock_info->renumber(i)] = indices_org[i]; + set_block_indices(); } } else @@ -79,8 +87,7 @@ namespace MeshWorker { indices_org.resize(c->get_fe().dofs_per_cell); c->get_mg_dof_indices(indices_org); - for (unsigned int i=0;iblock_info->renumber(i)] = indices_org[i]; + set_block_indices(); } level_cell = true; } diff --git a/deal.II/include/deal.II/meshworker/simple.h b/deal.II/include/deal.II/meshworker/simple.h index f3c5e6b4c2..6aec7dcd69 100644 --- a/deal.II/include/deal.II/meshworker/simple.h +++ b/deal.II/include/deal.II/meshworker/simple.h @@ -301,6 +301,38 @@ namespace MeshWorker */ void initialize(const MGConstrainedDoFs& mg_constrained_dofs); + /** + * Store information on the + * local block structure. If + * the assembler is + * inititialized with this + * function, + * initialize_info() will + * generate one local matrix + * for each block row and + * column, which will be + * numbered + * lexicographically, row by + * row. + * + * In spite of using local + * block structure, all + * blocks will be enteres + * into the same global + * matrix, disregarding any + * global block structure. + * + * @note The argument of this + * function will be copied + * into the member object + * #local_indices. Thus, + * every subsequent change in + * the block structure must + * be initialzied or will not + * be used by the assembler. + */ + void initialize_local_blocks(const BlockIndices& local_indices); + /** * Initialize the matrices * #flux_up and #flux_down @@ -461,6 +493,16 @@ namespace MeshWorker */ SmartPointer > mg_constrained_dofs; + /** + * The object containing the + * local block structure. Set + * by + * initialize_local_blocks() + * and used by assembling + * functions. + */ + BlockIndices local_indices; + /** * The smallest positive * number that will be @@ -768,6 +810,14 @@ namespace MeshWorker mg_constrained_dofs = &c; } + template + inline void + MGMatrixSimple::initialize_local_blocks(const BlockIndices& b) + { + local_indices = b; + } + + template inline void MGMatrixSimple::initialize_fluxes( @@ -793,7 +843,26 @@ namespace MeshWorker inline void MGMatrixSimple::initialize_info(DOFINFO& info, bool face) const { - info.initialize_matrices(1, face); + const unsigned int n = local_indices.size(); + + if (n == 0) + info.initialize_matrices(1, face); + else + { + info.initialize_matrices(n*n, face); + unsigned int k=0; + for (unsigned int i=0;i::assemble(const DOFINFO& info) { const unsigned int level = info.cell->level(); - assemble((*matrix)[level], info.matrix(0,false).matrix, info.indices, info.indices, level); - if(mg_constrained_dofs != 0) + if (local_indices.size() == 0) { - assemble_in((*interface_in)[level], info.matrix(0,false).matrix, info.indices, info.indices, level); - assemble_out((*interface_out)[level],info.matrix(0,false).matrix, info.indices, info.indices, level); + assemble((*matrix)[level], info.matrix(0,false).matrix, + info.indices, info.indices, level); + if(mg_constrained_dofs != 0) + { + assemble_in((*interface_in)[level], info.matrix(0,false).matrix, + info.indices, info.indices, level); + assemble_out((*interface_out)[level],info.matrix(0,false).matrix, + info.indices, info.indices, level); + } } + else + for (unsigned int k=0;k template inline void @@ -1044,39 +1136,79 @@ namespace MeshWorker const unsigned int level1 = info1.cell->level(); const unsigned int level2 = info2.cell->level(); - if (level1 == level2) + + if (local_indices.size() == 0) { - if(mg_constrained_dofs == 0) + if (level1 == level2) { - assemble((*matrix)[level1], info1.matrix(0,false).matrix, info1.indices, info1.indices); - assemble((*matrix)[level1], info1.matrix(0,true).matrix, info1.indices, info2.indices); - assemble((*matrix)[level1], info2.matrix(0,false).matrix, info2.indices, info2.indices); - assemble((*matrix)[level1], info2.matrix(0,true).matrix, info2.indices, info1.indices); + if(mg_constrained_dofs == 0) + { + assemble((*matrix)[level1], info1.matrix(0,false).matrix, info1.indices, info1.indices); + assemble((*matrix)[level1], info1.matrix(0,true).matrix, info1.indices, info2.indices); + assemble((*matrix)[level1], info2.matrix(0,false).matrix, info2.indices, info2.indices); + assemble((*matrix)[level1], info2.matrix(0,true).matrix, info2.indices, info1.indices); + } + else + { + assemble((*matrix)[level1], info1.matrix(0,false).matrix, info1.indices, info1.indices, level1); + assemble((*matrix)[level1], info1.matrix(0,true).matrix, info1.indices, info2.indices, level1); + assemble((*matrix)[level1], info2.matrix(0,false).matrix, info2.indices, info2.indices, level1); + assemble((*matrix)[level1], info2.matrix(0,true).matrix, info2.indices, info1.indices, level1); + } } - else + else { - assemble((*matrix)[level1], info1.matrix(0,false).matrix, info1.indices, info1.indices, level1); - assemble((*matrix)[level1], info1.matrix(0,true).matrix, info1.indices, info2.indices, level1); - assemble((*matrix)[level1], info2.matrix(0,false).matrix, info2.indices, info2.indices, level1); - assemble((*matrix)[level1], info2.matrix(0,true).matrix, info2.indices, info1.indices, level1); + Assert(level1 > level2, ExcInternalError()); + // Do not add info2.M1, + // which is done by + // the coarser cell + assemble((*matrix)[level1], info1.matrix(0,false).matrix, info1.indices, info1.indices); + if(level1>0) + { + assemble_up((*flux_up)[level1],info1.matrix(0,true).matrix, info2.indices, info1.indices, level1); + assemble_down((*flux_down)[level1], info2.matrix(0,true).matrix, info2.indices, info1.indices, level1); + } } } else + for (unsigned int k=0;k level2, ExcInternalError()); - // Do not add info2.M1, - // which is done by - // the coarser cell - assemble((*matrix)[level1], info1.matrix(0,false).matrix, info1.indices, info1.indices); - //assemble_transpose((*flux_up)[level1],info1.matrix(0,true).matrix, info2.indices, info1.indices); - //assemble((*flux_down)[level1], info2.matrix(0,true).matrix, info2.indices, info1.indices); - if(level1>0) + const unsigned int row = info1.matrix(k,false).row; + const unsigned int column = info1.matrix(k,false).column; + + if (level1 == level2) + { + if(mg_constrained_dofs == 0) + { + assemble((*matrix)[level1], info1.matrix(k,false).matrix, info1.indices_by_block[row], info1.indices_by_block[column]); + assemble((*matrix)[level1], info1.matrix(k,true).matrix, info1.indices_by_block[row], info2.indices_by_block[column]); + assemble((*matrix)[level1], info2.matrix(k,false).matrix, info2.indices_by_block[row], info2.indices_by_block[column]); + assemble((*matrix)[level1], info2.matrix(k,true).matrix, info2.indices_by_block[row], info1.indices_by_block[column]); + } + else + { + assemble((*matrix)[level1], info1.matrix(k,false).matrix, info1.indices_by_block[row], info1.indices_by_block[column], level1); + assemble((*matrix)[level1], info1.matrix(k,true).matrix, info1.indices_by_block[row], info2.indices_by_block[column], level1); + assemble((*matrix)[level1], info2.matrix(k,false).matrix, info2.indices_by_block[row], info2.indices_by_block[column], level1); + assemble((*matrix)[level1], info2.matrix(k,true).matrix, info2.indices_by_block[row], info1.indices_by_block[column], level1); + } + } + else { - assemble_up((*flux_up)[level1],info1.matrix(0,true).matrix, info2.indices, info1.indices, level1); - assemble_down((*flux_down)[level1], info2.matrix(0,true).matrix, info2.indices, info1.indices, level1); + Assert(level1 > level2, ExcInternalError()); + // Do not add info2.M1, + // which is done by + // the coarser cell + assemble((*matrix)[level1], info1.matrix(k,false).matrix, info1.indices_by_block[row], info1.indices_by_block[column]); + if(level1>0) + { + assemble_up((*flux_up)[level1],info1.matrix(k,true).matrix, info2.indices_by_block[row], info1.indices_by_block[column], level1); + assemble_down((*flux_down)[level1], info2.matrix(k,true).matrix, info2.indices_by_block[row], info1.indices_by_block[column], level1); + } } } } + //----------------------------------------------------------------------// template