From: bangerth Date: Tue, 4 Sep 2012 13:46:12 +0000 (+0000) Subject: Fix a bug where we try to access elements in a matrix that are beyond the end. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c0f317606bf5c3b48aacaeeb7626b6eaa010220e;p=dealii-svn.git Fix a bug where we try to access elements in a matrix that are beyond the end. git-svn-id: https://svn.dealii.org/trunk@26233 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 56a6582e34..62ef820261 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -50,7 +50,11 @@ never working correctly and it is not used.

Specific improvements

    -
  1. Nothing so far. +
  2. Fixed: When applying a ConstraintMatrix to a block matrix +where the last few rows are empty, we ran into an unrelated assertion. +This is now fixed. +
    +(Jason Sheldon, Wolfgang Bangerth, 2012/09/04)
diff --git a/deal.II/include/deal.II/base/thread_management.h b/deal.II/include/deal.II/base/thread_management.h index 6b5f941cfc..d2f83f0333 100644 --- a/deal.II/include/deal.II/base/thread_management.h +++ b/deal.II/include/deal.II/base/thread_management.h @@ -3739,7 +3739,6 @@ namespace Threads template struct TaskDescriptor; - /** * The task class for TBB that is * used by the TaskDescriptor @@ -3910,7 +3909,7 @@ namespace Threads template friend struct TaskEntryPoint; - template friend struct Task; + friend class dealii::Threads::Task; }; diff --git a/deal.II/include/deal.II/lac/constraint_matrix.templates.h b/deal.II/include/deal.II/lac/constraint_matrix.templates.h index ea3f98ab40..db31f44eb4 100644 --- a/deal.II/include/deal.II/lac/constraint_matrix.templates.h +++ b/deal.II/include/deal.II/lac/constraint_matrix.templates.h @@ -1768,9 +1768,9 @@ namespace internals const unsigned int loc_row = global_rows.local_row(i); #ifdef DEBUG - const unsigned int * col_ptr = sparsity.n_nonzero_elements() == 0 ? 0 : + const unsigned int * col_ptr = sparsity.row_length(row) == 0 ? 0 : &sparsity_struct[row_start[row]]; - number * val_ptr = sparsity.n_nonzero_elements() == 0 ? 0 : + number * val_ptr = sparsity.row_length(row) == 0 ? 0 : &sparse_matrix->global_entry (row_start[row]); #else const unsigned int * col_ptr = &sparsity_struct[row_start[row]]; diff --git a/tests/deal.II/constraints_block_01.cc b/tests/deal.II/constraints_block_01.cc new file mode 100644 index 0000000000..e2ed848947 --- /dev/null +++ b/tests/deal.II/constraints_block_01.cc @@ -0,0 +1,348 @@ +//---------------------------- constraints.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2012 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- constraints.cc --------------------------- + + +// test that we fixed a case where we tried to deal with some constraints in a +// block matrix where the last few rows of one of the blocks were empty and we +// ran into an unrelated assertion because we were accessing something beyond +// the end of the array +// +// testcase by Jason Sheldon + +#include "../tests.h" +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include + + +std::ofstream logfile("constraints_block_01/output"); + + +int main () +{ + deallog << std::setprecision (2); + logfile << std::setprecision (2); + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + // setting up some constants + const unsigned int dim = 2; + const unsigned int solid_dim = 2*dim; + const unsigned int fluid_dim = dim+1; + const unsigned int mesh_dim = dim; + const unsigned int total_dim = solid_dim + fluid_dim + mesh_dim; + + // make the tria and domain + Triangulation tria; + + GridGenerator::hyper_cube (tria); + tria.refine_global(1); + + Point real_cell_center; + + // create a solid domain (0) in the + // lower left hand corner of a 2x2 grid + // the rest is fluid/mesh (1) + + for (Triangulation::active_cell_iterator + cell = tria.begin_active(); + cell != tria.end(); ++cell) + { + real_cell_center = cell->center(); + if (real_cell_center(0) < 0.5 + && real_cell_center(1) < 0.5 ) + cell->set_material_id (0); // solid + else + cell->set_material_id (1); //fluid + }//cell + + //create the FE spaces for the solid and the fluid/mesh + //each are padded with FE_Nothing to be equal length + + std::string solid_fe_name = "FESystem[FE_Q(2)^2-FE_Q(2)^2-FE_Nothing()^2-FE_Nothing()-FE_Nothing()^2]"; + std::string fluid_fe_name = "FESystem[FE_Nothing()^2-FE_Nothing()^2-FE_Q(2)^2-FE_Q(1)-FE_Q(2)^2]"; + + hp::FECollection fe_collection; + FiniteElement * solid_fe = FETools::get_fe_from_name(solid_fe_name); + FiniteElement * fluid_fe = FETools::get_fe_from_name(fluid_fe_name); + + deallog << "Solid FE Space: " << solid_fe->get_name() << std::endl; + deallog << "Fluid/Mesh FE Space: " << fluid_fe->get_name() << std::endl; + + fe_collection.push_back(*solid_fe); + fe_collection.push_back(*fluid_fe); + + hp::DoFHandler dh(tria); + + for (hp::DoFHandler::active_cell_iterator + cell = dh.begin_active(); + cell != dh.end(); ++cell) + { + if (int(cell->material_id()) == 0) + cell->set_active_fe_index (0); + else + cell->set_active_fe_index (1); + }//cell + dh.distribute_dofs(fe_collection); + + std::vector block_component (total_dim,0); + + for (unsigned int comp=0; comp < total_dim; comp++) + { + if ( comp < solid_dim ) + block_component[comp] = 0; + else if ( comp < solid_dim+fluid_dim ) + block_component[comp] = 1; + else + block_component[comp] = 2; + }//comp + + std::vector dofs_per_block(3,0);//3 blocks, count dofs: + DoFTools::count_dofs_per_component (dh, dofs_per_block, false, block_component); + + DoFRenumbering::component_wise(dh, block_component); + + //build the sparsitypattern + + BlockSparsityPattern block_sparsity_pattern; + { + BlockCompressedSimpleSparsityPattern csp(3,3); + csp.block(0,0).reinit (dofs_per_block[0],dofs_per_block[0]);//solid-solid + csp.block(0,1).reinit (dofs_per_block[0],dofs_per_block[1]);//solid-fluid + csp.block(0,2).reinit (dofs_per_block[0],dofs_per_block[2]);//solid-mesh + + csp.block(1,0).reinit (dofs_per_block[1],dofs_per_block[0]);//fluid-solid + csp.block(1,1).reinit (dofs_per_block[1],dofs_per_block[1]);//fluid-fluid + csp.block(1,2).reinit (dofs_per_block[1],dofs_per_block[2]);//fluid-mesh + + csp.block(2,0).reinit (dofs_per_block[2],dofs_per_block[0]);//mesh-solid + csp.block(2,1).reinit (dofs_per_block[2],dofs_per_block[1]);//mesh-fluid + csp.block(2,2).reinit (dofs_per_block[2],dofs_per_block[2]);//mesh-mesh + + csp.collect_sizes(); + + //enforce coupling across cells and interface + + Table<2,DoFTools::Coupling> cell_coupling (fe_collection.n_components(), fe_collection.n_components()); + Table<2,DoFTools::Coupling> face_coupling (fe_collection.n_components(), fe_collection.n_components()); + + for (unsigned int c=0; c=solid_dim) && (d>=solid_dim)) //couples fluid dims and mesh dims + && !((c==solid_dim+dim) && (d==solid_dim+dim)))) //fluid pressure does not couple with itself + { + cell_coupling[c][d] = DoFTools::always; + cell_coupling[d][c] = DoFTools::always; + }//cell_coupling + + if ((c=solid_dim)) //couples solid dims with fluid and mesh dims on the interface + { + face_coupling[c][d] = DoFTools::always; + face_coupling[d][c] = DoFTools::always; + }//face_coupling + }//d + }//c + DoFTools::make_flux_sparsity_pattern (dh, csp, cell_coupling, face_coupling); + block_sparsity_pattern.copy_from (csp); + } + + // build matrices and vectors + + BlockSparseMatrix system_matrix(block_sparsity_pattern); + BlockVector system_rhs(dofs_per_block); + + // set up constraints for solution and solution update + /** + * These constraints enforce that the solid displacement + * and mesh displacement are the same on the interface. + * + * They also constrain that the solid velocity + * and fluid velocity are the same on the interface + * + * The interface check is simplified for this 2x2 case + */ + + ConstraintMatrix constraints; + + const unsigned int dofs_per_fl_msh_face = fluid_fe->dofs_per_face; + const unsigned int dofs_per_solid_face = solid_fe->dofs_per_face; + std::vector fl_msh_face_dof_indices (dofs_per_fl_msh_face); + std::vector solid_face_dof_indices (dofs_per_solid_face ); + + std::vector > solid_fluid_pairs; + std::vector > solid_mesh_pairs; + + for (hp::DoFHandler::active_cell_iterator + cell = dh.begin_active(); cell!=dh.end(); ++cell) //loops over the cells + { + if (int(cell->material_id()) == 1) + {//Only loop over cells in the fluid region + for (unsigned int f=0; f::faces_per_cell; ++f) + { + if (!cell->at_boundary(f)) + { + bool face_is_on_interface = false; + //checks to see if cell face neighbor is in the solid domain + if (int(cell->neighbor(f)->material_id()) == 0 ) + face_is_on_interface = true; + + if (face_is_on_interface) + { + std::vector solid_disp_dof, solid_vel_dof, fluid_vel_dof, mesh_disp_dof; + + cell->face(f)->get_dof_indices (solid_face_dof_indices, 0); + cell->face(f)->get_dof_indices (fl_msh_face_dof_indices, 1); + + for (unsigned int i=0; iface_system_to_component_index(i).first; + if (comp = dim && comp < solid_dim) + solid_vel_dof.push_back(i); + }//i + for (unsigned int i=0; iface_system_to_component_index(i).first; + if (comp >= solid_dim && comp < solid_dim+dim) + fluid_vel_dof.push_back(i); + if (comp>=solid_dim+fluid_dim) + mesh_disp_dof.push_back(i); + }//i + + for (unsigned int i=0; i + (solid_face_dof_indices[solid_vel_dof[i]], + fl_msh_face_dof_indices[fluid_vel_dof[i]])); + }//i + for (unsigned int i=0; i + (solid_face_dof_indices[solid_disp_dof[i]], + fl_msh_face_dof_indices[mesh_disp_dof[i]])); + } + }//at interface check + }//not at boundary check + }//face + }//is this in the fluid material region? + }//cell + + constraints.close(); + + // prints out which dofs are coupled + deallog<<"---------------Coupled dofs---------------"<dofs_per_cell; + + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + Vector cell_rhs (dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + for (hp::DoFHandler::active_cell_iterator + cell = dh.begin_active(); + cell!=dh.end(); + ++cell) //loops over the cells + { + if (int(cell->material_id()) == 1)//are we in the fluid region? + { + cell->get_dof_indices (local_dof_indices); + + // normally the governing equations would be added + // to the cell matrix and right hand side here + + deallog<<"I am about to distribute cell: "<center()<[FE_Q<2>(2)^2-FE_Q<2>(2)^2-FE_Nothing<2>()^2-FE_Nothing<2>()-FE_Nothing<2>()^2] +DEAL::Fluid/Mesh FE Space: FESystem<2>[FE_Nothing<2>()^2-FE_Nothing<2>()^2-FE_Q<2>(2)^2-FE_Q<2>(1)-FE_Q<2>(2)^2] +DEAL::---------------Coupled dofs--------------- +DEAL::solid dof: 6, fluid dof: 36 +DEAL::solid dof: 7, fluid dof: 37 +DEAL::solid dof: 14, fluid dof: 42 +DEAL::solid dof: 15, fluid dof: 43 +DEAL::solid dof: 22, fluid dof: 48 +DEAL::solid dof: 23, fluid dof: 49 +DEAL::solid dof: 10, fluid dof: 58 +DEAL::solid dof: 11, fluid dof: 59 +DEAL::solid dof: 14, fluid dof: 42 +DEAL::solid dof: 15, fluid dof: 43 +DEAL::solid dof: 30, fluid dof: 71 +DEAL::solid dof: 31, fluid dof: 72 +DEAL::solid dof: 4, mesh dof: 86 +DEAL::solid dof: 5, mesh dof: 87 +DEAL::solid dof: 12, mesh dof: 90 +DEAL::solid dof: 13, mesh dof: 91 +DEAL::solid dof: 20, mesh dof: 94 +DEAL::solid dof: 21, mesh dof: 95 +DEAL::solid dof: 8, mesh dof: 104 +DEAL::solid dof: 9, mesh dof: 105 +DEAL::solid dof: 12, mesh dof: 90 +DEAL::solid dof: 13, mesh dof: 91 +DEAL::solid dof: 28, mesh dof: 114 +DEAL::solid dof: 29, mesh dof: 115 +DEAL::------------------------------------------ +DEAL::I am about to distribute cell: 1.1 +DEAL::Cell: 1.1 has it's center at 0.75 0.25 +DEAL::I am about to distribute cell: 1.2 +DEAL::Cell: 1.2 has it's center at 0.25 0.75 +DEAL::I am about to distribute cell: 1.3 +DEAL::Cell: 1.3 has it's center at 0.75 0.75