From: Baerbel Jannsen Date: Tue, 14 Sep 2010 14:53:00 +0000 (+0000) Subject: do some modifications to be able to use multigrid for Raviart Thomas X-Git-Tag: v8.0.0~5515 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=86241fcce20262d405e30bd124f4579644ce84bf;p=dealii.git do some modifications to be able to use multigrid for Raviart Thomas elements git-svn-id: https://svn.dealii.org/trunk@21958 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/multigrid/mg_constrained_dofs.h b/deal.II/deal.II/include/multigrid/mg_constrained_dofs.h index e5d0045dce..9130535722 100644 --- a/deal.II/deal.II/include/multigrid/mg_constrained_dofs.h +++ b/deal.II/deal.II/include/multigrid/mg_constrained_dofs.h @@ -119,6 +119,12 @@ class MGConstrainedDoFs : public Subscriptor const std::vector > & get_refinement_edge_boundary_indices () const; + /** + * Return if boundary_indices need to + * be set or not. + */ + + const bool set_boundary_values () const; private: /** @@ -263,6 +269,16 @@ MGConstrainedDoFs::get_refinement_edge_boundary_indices () const return refinement_edge_boundary_indices; } +inline +const bool +MGConstrainedDoFs::set_boundary_values () const +{ + const bool boundary_values_need_to_be_set + = boundary_indices.size()!=0; + return boundary_values_need_to_be_set; +} + + DEAL_II_NAMESPACE_CLOSE #endif diff --git a/deal.II/deal.II/include/numerics/mesh_worker_assembler.h b/deal.II/deal.II/include/numerics/mesh_worker_assembler.h index dddff77ecb..1c29f66813 100644 --- a/deal.II/deal.II/include/numerics/mesh_worker_assembler.h +++ b/deal.II/deal.II/include/numerics/mesh_worker_assembler.h @@ -574,14 +574,14 @@ namespace MeshWorker /** * Initialize the matrices - * #interface_up and #interface_down + * #interface_in and #interface_out * used for local refinement * with continuous * Galerkin methods. */ - void initialize_interfaces(MGLevelObject& interface_up, - MGLevelObject& interface_down); + void initialize_interfaces(MGLevelObject& interface_in, + MGLevelObject& interface_out); /** * Initialize the local data * in the @@ -626,10 +626,19 @@ namespace MeshWorker * into a global matrix. */ void assemble(MATRIX& G, + const FullMatrix& M, + const std::vector& i1, + const std::vector& i2); + + /** + * Assemble a single matrix + * into a global matrix. + */ + void assemble(MATRIX& G, const FullMatrix& M, const std::vector& i1, const std::vector& i2, - const unsigned int level = -1); + const unsigned int level); /** * Assemble a single matrix @@ -638,14 +647,23 @@ namespace MeshWorker void assemble_transpose(MATRIX& G, const FullMatrix& M, const std::vector& i1, - const std::vector& i2, - const unsigned int level = -1); + const std::vector& i2); /** * Assemble a single matrix * into a global matrix. */ - void assemble_down(MATRIX& G, + void assemble_in(MATRIX& G, + const FullMatrix& M, + const std::vector& i1, + const std::vector& i2, + const unsigned int level = -1); + + /** + * Assemble a single matrix + * into a global matrix. + */ + void assemble_out(MATRIX& G, const FullMatrix& M, const std::vector& i1, const std::vector& i2, @@ -680,7 +698,7 @@ namespace MeshWorker * refinement edge, coupling * coarse to fine. */ - SmartPointer,MGMatrixSimple > interface_up; + SmartPointer,MGMatrixSimple > interface_in; /** * The matrix used for face @@ -689,7 +707,7 @@ namespace MeshWorker * refinement edge, coupling * fine to coarse. */ - SmartPointer,MGMatrixSimple > interface_down; + SmartPointer,MGMatrixSimple > interface_out; /** * A pointer to the object containing constraints. */ @@ -1463,10 +1481,10 @@ namespace MeshWorker template inline void MGMatrixSimple::initialize_interfaces( - MGLevelObject& up, MGLevelObject& down) + MGLevelObject& in, MGLevelObject& out) { - interface_up = &up; - interface_down = &down; + interface_in = ∈ + interface_out = &out; } @@ -1480,6 +1498,23 @@ namespace MeshWorker } + template + inline void + MGMatrixSimple::assemble( + MATRIX& G, + const FullMatrix& M, + const std::vector& i1, + const std::vector& i2) + { + AssertDimension(M.m(), i1.size()); + AssertDimension(M.n(), i2.size()); + + for (unsigned int j=0; j= threshold) + G.add(i1[j], i2[k], M(j,k)); + } + template inline void @@ -1505,23 +1540,46 @@ namespace MeshWorker for (unsigned int j=0; j= threshold) - if(mg_constrained_dofs->at_refinement_edge(level, i1[j])==false && - mg_constrained_dofs->at_refinement_edge(level, i2[k])==false) - if((mg_constrained_dofs->is_boundary_index(level, i1[j])==false && - mg_constrained_dofs->is_boundary_index(level, i2[k])==false) - || - (mg_constrained_dofs->is_boundary_index(level, i1[j])==true && - mg_constrained_dofs->is_boundary_index(level, i2[k])==true - && - i1[j] == i2[k]) - ) - G.add(i1[j], i2[k], M(j,k)); + if (!mg_constrained_dofs->at_refinement_edge(level, i1[j]) && + !mg_constrained_dofs->at_refinement_edge(level, i2[k])) + { + if (mg_constrained_dofs->set_boundary_values()) + { + if ((!mg_constrained_dofs->is_boundary_index(level, i1[j]) && + !mg_constrained_dofs->is_boundary_index(level, i2[k])) + || + (mg_constrained_dofs->is_boundary_index(level, i1[j]) && + mg_constrained_dofs->is_boundary_index(level, i2[k]) && + i1[j] == i2[k])) + G.add(i1[j], i2[k], M(j,k)); + } + else + G.add(i1[j], i2[k], M(j,k)); + } } } + + template + inline void + MGMatrixSimple::assemble_transpose( + MATRIX& G, + const FullMatrix& M, + const std::vector& i1, + const std::vector& i2) + { + AssertDimension(M.n(), i1.size()); + AssertDimension(M.m(), i2.size()); + + for (unsigned int j=0; j= threshold) + G.add(i1[j], i2[k], M(k,j)); + } + template inline void - MGMatrixSimple::assemble_down( + MGMatrixSimple::assemble_in( MATRIX& G, const FullMatrix& M, const std::vector& i1, @@ -1543,23 +1601,28 @@ namespace MeshWorker for (unsigned int j=0; j= threshold) - if(mg_constrained_dofs->at_refinement_edge(level, i1[j])==true && - mg_constrained_dofs->at_refinement_edge(level, i2[k])==false) - if((mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j])==false && - mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k])==false) - || - (mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j])==true && - mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k])==true - && - i1[j] == i2[k]) - ) - G.add(i1[j], i2[k], M(j,k)); + if(mg_constrained_dofs->at_refinement_edge(level, i1[j]) && + !mg_constrained_dofs->at_refinement_edge(level, i2[k])) + { + if (mg_constrained_dofs->set_boundary_values()) + { + if((!mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) && + !mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k])) + || + (mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) && + mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]) && + i1[j] == i2[k])) + G.add(i1[j], i2[k], M(j,k)); + } + else + G.add(i1[j], i2[k], M(j,k)); + } } } template inline void - MGMatrixSimple::assemble_transpose( + MGMatrixSimple::assemble_out( MATRIX& G, const FullMatrix& M, const std::vector& i1, @@ -1581,17 +1644,22 @@ namespace MeshWorker for (unsigned int j=0; j= threshold) - if(mg_constrained_dofs->at_refinement_edge(level, i1[j])==true && - mg_constrained_dofs->at_refinement_edge(level, i2[k])==false) - if((mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j])==false && - mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k])==false) - || - (mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j])==true && - mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k])==true - && - i1[j] == i2[k]) - ) - G.add(i1[j], i2[k], M(k,j)); + if(mg_constrained_dofs->at_refinement_edge(level, i1[j]) && + !mg_constrained_dofs->at_refinement_edge(level, i2[k])) + { + if (mg_constrained_dofs->set_boundary_values()) + { + if((!mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) && + !mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k])) + || + (mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) && + mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]) && + i1[j] == i2[k])) + G.add(i1[j], i2[k], M(k,j)); + } + else + G.add(i1[j], i2[k], M(k,j)); + } } } @@ -1605,8 +1673,8 @@ namespace MeshWorker assemble((*matrix)[level], info.matrix(0,false).matrix, info.indices, info.indices, level); if(mg_constrained_dofs != 0) { - assemble_transpose((*interface_down)[level],info.matrix(0,false).matrix, info.indices, info.indices, level); - assemble_down((*interface_up)[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_in((*interface_in)[level], info.matrix(0,false).matrix, info.indices, info.indices, level); } } @@ -1629,6 +1697,13 @@ namespace MeshWorker 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 { diff --git a/deal.II/deal.II/source/multigrid/mg_transfer_prebuilt.cc b/deal.II/deal.II/source/multigrid/mg_transfer_prebuilt.cc index 6a02a3d6db..d0f13cfae5 100644 --- a/deal.II/deal.II/source/multigrid/mg_transfer_prebuilt.cc +++ b/deal.II/deal.II/source/multigrid/mg_transfer_prebuilt.cc @@ -285,8 +285,8 @@ void MGTransferPrebuilt::build_matrices ( // impose boundary conditions // but only in the column of // the prolongation matrix - if(mg_constrained_dofs != 0) - if (mg_constrained_dofs->get_boundary_indices().size() != 0) + if (mg_constrained_dofs != 0) + if (mg_constrained_dofs->set_boundary_values()) { std::vector constrain_indices; for (int level=n_levels-2; level>=0; --level)